Computational complexity continuum within Ising formulation of NP problems

A promising approach to achieve computational supremacy over the classical von Neumann architecture explores classical and quantum hardware as Ising machines. The minimisation of the Ising Hamiltonian is known to be NP-hard problem yet not all problem instances are equivalently hard to optimise. Given that the operational principles of Ising machines are suited to the structure of some problems but not others, we propose to identify computationally simple instances with an ‘optimisation simplicity criterion’. Neuromorphic architectures based on optical, photonic, and electronic systems can naturally operate to optimise instances satisfying this criterion, which are therefore often chosen to illustrate the computational advantages of new Ising machines. As an example, we show that the Ising model on the Möbius ladder graph is ‘easy’ for Ising machines. By rewiring the Möbius ladder graph to random 3-regular graphs, we probe an intermediate computational complexity between P and NP-hard classes with several numerical methods. Significant fractions of polynomially simple instances are further found for a wide range of small size models from spin glasses to maximum cut problems. A compelling approach for distinguishing easy and hard instances within the same NP-hard class of problems can be a starting point in developing a standardised procedure for the performance evaluation of emerging physical simulators and physics-inspired algorithms. The advantage of unconventional computing architectures is commonly demonstrated by solving an NP-hard problem, but some instances are easier to solve than others. Here, an optimisation simplicity criterion is proposed that classifies the complexity of instances on optical or electronic neuromorphic computers.

T he recent advances in developing physical platforms for solving combinatorial optimisation problems reveal the future of high-performance computing for quantum and classical devices.Unconventional computing architectures were proposed for numerous systems, including superconducting qubits [1][2][3] , field-programmable gate arrays 4 , optical parametric oscillators 5,6 , memristors 7 , lasers [8][9][10] , photonic simulators 11,12 , trapped ions 13 , polariton 14,15 and photon 16 condensates.An attractive opportunity to show the advantageous performance of one system over others becomes a demonstration of the platform's ability to optimise non-deterministic polynomial time (NP-hard) problems that are computationally intractable for the traditional von Neumann architecture machines.The intractability is manifested in the exponential growth of the number of operations with the problem size.From the computational complexity theory perspective, the exponential growth does not necessarily apply to all instances of an optimisation problem, which is shown to be NP-hard in general, admitting the worstcase scenario when a mere handful of instances are truly hard to optimise.Selection of the hardest instances within NP-hard classes could be the key to determining the computational advantages of small and medium-size simulators and may lead to a reliable generalisation of their optimisation performance to a larger scale.
The hard optimisation problems from vastly different areas, including the travelling salesman problem, spin glass models, knapsack problem, integer linear programming, can be reformulated to minimise spin Hamiltonians 17 , among which a special place belongs to the Ising Hamiltonian.To minimise the Ising Hamiltonian ('solve the Ising model'), one needs to find N binary spins s i ∈ {−1, 1} that minimise where J ij are real coupling coefficients and h i are external fields.
Solving the Ising model is NP-hard problem in general, with computational hardness proven for certain coupling matrices 18 (see Methods for details).The Ising Hamiltonian is universal, meaning a fine-graining procedure exists to transform any classical spin Hamiltonian, continuous or discrete, with an arbitrary coupling matrix to the low energy spectra of the universal model such as the Ising model on the square lattice with fields 19 .
Given existing small and medium-scale simulators, considerable attention is devoted to problems that can be mapped to the Ising model with zero overhead.A common example includes the maximum cut (MaxCut) class of problems in which one looks for the cut of the given graph into two subsets with the largest number of their connecting weighted edges.The subclass of unweighted graphs is attractive for experimental implementation since it only requires the realisation of antiferromagnetic couplings (J ij < 0) of the same amplitude, i.e., J ij = −1 if spins i and j are connected, and 0 otherwise.Since the unweighted MaxCut problem is NP-hard 20 , the instances of unweighted k-regular graphs, in which each spin is randomly connected to k other spins, are often used to study new and compare existing physical simulators 5,11,12,21,22 .The 3-regular MaxCut problems are used in the proposal of the quantum approximate optimisation algorithm 23 with its later experimental demonstration on superconducting qubits 24 .
Another common practice is to consider the unweighted MaxCut problems on circulant graphs.Circulant graphs are defined by symmetric circulant adjacency matrices where (i + 1)th row is a cyclic shift of i-th row by one element to the right.Subclasses of circulant graphs include complete graphs, cyclic graphs, Möbius ladder and many others 25,26 .Efficient quantum walks are implemented on circulant graphs with sampling problems shown to be intractable for classical hardware 27 .The complete unweighted graphs with antiferromagnetic couplings can be optimised for large sizes up to 80,000 with the photonic Ising machine 11 .The unweighted Ising model on the Möbius ladder graph formally belongs to the MaxCut problem and has the circulant adjacency matrix with nonzero elements of the first row at 0, N/2, and N-th positions, where N is an even number.For the Möbius ladder of size N = 100, the ground state can be found with a probability of 21% for the coherent Ising machine based on optical parametric oscillators 5,28,29 and a success rate of 34% is demonstrated with optoelectronic oscillators 30 .The 3% ground state probability is reported for the larger Möbius ladder of size 300 on the analogue coupled electronic oscillator machine 31 .The Möbius ladders become typical candidates for evaluating the performance of physical platforms 32,33 and an exponential time increase on the graphs up to 800 nodes has been reported 34 .Ordinarily, it is tempting to assume that choosing any instance of a general class of NP-hard problems is equivalent to considering a hard instance, thereby ignoring the possibility of that instance being in the P-class.
In this article, we probe an instance complexity between the two extremes.To detect easy instances within the Ising model, we propose an 'optimisation simplicity criterion'.We provide numerical evidence of such optimisation simplicity for instances covering a wide range of problems from spin glass models to 3-regular MaxCut problems.As an illustrative example of easy instances of the unweighted 3-regular MaxCut problem, the Möbius ladder graphs are shown to be polynomially solvable.In particular, greater than 99% ground state probability can be ensured with the quadratic increase in the number of time iterations for the HT algorithm 35 on graphs up to 10,000 size.Moreover, the mathematical complexity of the weighted Ising model on the Möbius ladder graphs is shown to be in P-class, and the superlinear scaling for its computational complexity is demonstrated with the exact commercial solver, Gurobi.With a simple Möbius ladder at one end and hard random 3-regular graph at the other, the relative computational hardness of intermediate graphs with rearranged edges is investigated.The percentage of rewired edges in the Möbius ladder, which is required to achieve an average hardness of an arbitrary 3-regular graph, depends on the optimisation technique and can vary from 2-5% to 40-50%, as evidenced by the time performance of several heuristic algorithms and Gurobi solver.The Ising models satisfying the proposed optimisation simplicity criterion are not limited to circulant matrices and include sparse and dense interaction matrices of various topologies with or without a magnetic field.For some Ising models, such as the Mattis model, unweighted spin glasses on a torus, and biased ferromagnet on the Chimera graph, we find that all instances are solvable in polynomial time.There also exists a high probability of finding simple small size random instances of NP-hard problems, as we confirm for 3-regular MaxCut, Sherrington-Kirkpatrick, and other spin glass models, with couplings taken from the Gaussian and bimodal distributions.Understanding the average instance complexity of NP-hard problems and having a robust way to identify the polynomially easy instances could help evaluate the general potential of small and medium-scale simulators in solving hard combinatorial optimisation problems.

Results
The original work of Hopfield and Tank 35 introduced an analogue computational network for solving difficult optimisation problems.The network, later termed the Hopfield-Tank (HT) model or HT neural network, is governed by the equations: where x i (t) is a real input that describes the state of the i-th network element at time t, τ is the decay parameter, J is the symmetric coupling matrix, I b i are the offset biases (external fields) that can be absorbed into J by introducing an additional spin, N is the size of the network, and g(x i ) is the activation function.The nondecreasing monotonic function g(x i ) is designed to limit possible values of v i to the [−1, 1] range and is typically chosen as a sigmoid or hyperbolic tangent.The steady states of the HT model ( 2) are the minima of the Lyapunov function E: In the high-gain limit, when τ → ∞ or g approaches a step function g(x) = 1 (g(x) = −1) if x ≥ 0 (x < 0), the minima of E occur at v i = {−1, 1} and correspond to the minima of Eq. ( 1).If the high-gain limit conditions are violated (low-gain limit), the minima of E are not necessarily at v i = {−1, 1} and can be inside the hypercube [−1, 1] N .By projecting non-integer amplitudes of the steady-state at the end of the simulation, the allowed minimiser of the Ising model is restored at the nearest hypercube corner.Therefore, the HT network tends to locate local minima if it minimises the Ising model at all, as has been recognised in earlier works 36 .There exist simple coupling matrices J that can be globally optimised even in this low-gain limit.For zero fields in both limits, the steady states are completely characterised by the coupling matrix eigenvalues λ i and corresponding orthogonal eigenvectors e i 2 R N 1 with the coupling matrix expressed as In the presence of degenerate or zero eigenvalues, the eigenvectors form a subspace of rank lower than N. Denoting components of v in the space of coupling matrix eigenvectors as γ i and the null subspace component as q, the amplitudes and energy can be written as To minimise E, the components γ i should be increased for positive λ i and decreased otherwise.This observation reveals the nature of how the HT algorithm functions: it changes amplitudes v i in a way that gradually favours the larger positive eigenvalues λ i 37 .Therefore, in the low-gain limit, the HT algorithm finds the minimum of the Ising model that corresponds to the largest positive (dominant) eigenvalue.Suppose this minimum happens to be the global minimum, which is valid for many problems selected for testing the Ising Hamiltonian minimisers.In that case, the corresponding instances should be considered polynomially simple for optimisation, as we further explain in the paper.The choice of the HT algorithm in our analysis is not accidental and is prompted by its ability to replicate the dynamic behaviour of many existing Ising simulators considered in optics, photonics, and electronics.For instance, the recent memristorbased annealing system operates as a Hopfield neural network 7 .Another example is the coherent Ising machine on the optical parametric oscillators that is commonly thought to be similar to HT networks with nonlinear saturation of amplitudes and, therefore, both are often compared 21 .For such gain-dissipative computing machines, the successive better minima toward the dominant eigenmode 38 are achieved via a series of bifurcations 39 .
In general, the global minimum of the Ising Hamiltonian would correspond to a nontrivial direction in the eigenspace of e i in Eq. ( 5).This obvious yet substantial observation leads to our proposal for 'optimisation simplicity criterion': Proposition (Optimisation Simplicity Criterion-OSC).The instance of a hard problem should be regarded as computationally simple for optical and electronic analogue machines if the ground state minimiser s gs of the Ising Hamiltonian H Ising is located at the hypercube corner of the projected dominant eigenvector e max , corresponding to one of the largest eigenvalues λ max of the coupling matrix J: Without the loss of generality, the fields (the biases in HT networks) are assumed to be zero since they can always be incorporated into the coupling matrix J with an additional spin.The OSC provides an upper bound for the ground state energy of the Ising model.Eigenvalue analysis is common to many polynomial time algorithms that approximate both lower and upper bounds for optimal solutions to complex combinatorial problems 40,41 .For the MaxCut problem, the eigenvalue minimisation is known to be equivalent to semidefinite programming 42 , which in turn makes it equivalent to the eigenvalue maximisation that the HT algorithm does for the Ising model.
The standard procedure for verifying whether a particular instance satisfies the OSC would be to compare the upper bound energy E λ , which corresponds to the dominant eigenvector, with the global minimum obtained with a physical simulator or an optimisation algorithm.If these two energies coincide, the instance should be considered trivial to optimise.The polynomial complexity of instances satisfying the OSC could be recovered with the HT algorithm Eq. ( 2), which is naturally designed to project the input vector into a subspace that is dictated by the eigenvalues of the coupling matrix.For instance, to violate the OSC, it is sufficient to show that it has energy lower than E λ .The complexity of the instances that do not satisfy the OSC can be further assessed by other means.For example, the optimality gaps could be evaluated using the exact solvers such as Gurobi, or the time to solution metric could be considered for heuristic solvers, as we show below.
Möbius ladder graphs.As an illustrative example, we apply the HT algorithm Eq. ( 2) with a hyperbolic tangent activation function to a particular topology of unweighted 3-regular graphs, namely the Möbius ladder graph.The two representations of this cubic circulant graph of size N are shown in Fig. 1a.When n = N/2 is an even number, antiferromagnetic interactions cause lattice frustrations that result in N degenerate ground states with two frustrated edges (shown in red) between two domains of n anti-aligned spins and the ground state energy of (3n − 4). Figure 1b demonstrates a typical simulation of the HT network for the Möbius ladder of size N = 1000.The ratio of the HT energy, found by associating spins with the signs of amplitudes v i at the steady-state to the ground state energy is defined as the proximity to the ground state.The network operates in the low-gain limit (see Methods for parameters) and, hence, the amplitudes v i are not binary when the steady-state is reached.Yet, by gradually favouring the eigenvectors with larger eigenvalues, the HT algorithm moves spin states through the hypercube interior over time and achieves the global minimum, although the coupling matrix is modified by non-equal continuous amplitudes v i in [−1, 1].The necessity of homogeneous amplitudes for minimising nontrivial spin Hamiltonians with gaindissipative networks was discussed recently 43,44 .All states of the low energy spectra E λ i in Fig. 1b correspond to the eigenvectors of the largest eigenvalues of the interaction matrix, whose analytical expressions are available for the Möbius ladder as a representative of circulant matrices.
To estimate the performance of optical and electronic Ising machines, the Möbius ladder graphs, we determine the number of HT time iterations for achieving the ground state with probabilities greater than 50%, 70% and 99% for problem sizes up to N = 10,000.The ground state probability is defined as the fraction of simulations leading to the global minimum to the total number of simulations.Figure 1c shows a polynomial (quadratic) increase in the number of iterations with the graph size, which confirms the optimisation easiness of such problems.The quadratic slope remains the same for each range of the desired ground state probabilities.The ground state probability decreases for the fixed number of iterations, as demonstrated in Fig. 1d, which suggests that the reported quick performance deterioration of the physical Ising machines with the network size 5,30 may be caused by the fixed amount of internal system loops available in that physical platform.The lack of frustration in the Möbius ladders with odd N/2 does not necessarily mean that the ground state is trivial to reach.We consistently observe that such nonfrustrated graphs require larger number of time iterations than frustrated Möbius ladder graphs with even N/2.Since the complexity of one-time iteration of the HT algorithm is determined by the matrix-vector multiplication product as OðkNÞ for k-regular sparse graphs, the time complexity for globally optimising Ising Hamiltonian on the Möbius ladders scales as OðN 3 Þ with the problem size.
Thus, the eigenvalue maximisation principle, which underlies HT algorithm, ignores the energy profile of a simple problem that satisfies the OSC.Even in the absence of a mechanism for exploring the global energy landscape, network elements follow eigenvectors with successively larger eigenvalues.The corresponding consecutive energy states can differ by hundreds of spins, while the Ising Hamiltonian energy monotonically approaches the global minimum.This dynamic behaviour is drastically different from both the Ising machines based on equilibrium systems and optimisation methods, for which the width and height of energy barriers are critical and occasional increases in energy are common once the system escapes local minima.For example, the exponential time scaling for the Ising model on the Möbius ladders was recently reported for the simulated annealing algorithm 34 , while unconventional computing platforms based on gain-dissipative networks can efficiently find principal eigenvectors for million size problems 45 .
So far, we have discussed only the computational complexity of the unweighted Ising model on the Möbius ladder graphs.It can be seen that the mathematical complexity of minimising the Ising Hamiltonian on the Möbius ladder topology is in P-class.Since the Möbius ladder graph becomes a bipartite graph after the removal of two nodes, it belongs to a family of weakly-bipartite graphs 46 .However, the weighted MaxCut problem is in P-class for weakly-bipartite graphs 47 and, hence, the Ising model with arbitrary couplings on the Möbius ladder graph is in P-class too.
With the understanding of what is essential for an individual instance of the NP-hard problem to be counted as simple, we next present a natural approach for restoring complexity and study the continuous complexity transition from simple to hard instances for Ising optimisation on physical Ising machines.
Interpolating between simple Möbius ladders and hard 3-regular MaxCut.We develop a procedure that allows us to continuously 'tune' the graph from the Möbius ladders to random 3-regular graphs, the unweighted MaxCut problem on which is known to be NP-hard, and thereby to probe the intermediate problem computational complexity.To interpolate between two extremes, we consider the following random rewiring procedure.Starting from the Möbius ladder, we remove and reconnect a pair of edges at random.For each subsequent iteration of the rewiring procedure, a random pair among the original edges (if any) of the Möbius ladder is selected.Hence, intermediate graphs are quantified by the percentage of rewired edges in the Möbius ladder.For the frustrated Möbius ladder graphs to violate the OSC, the rearrangement of two edges is sufficient for any problem size N as shown in Fig. 2a and works for about 85% of the Möbius  ladders of size up to 1000 in Fig. 2b.Both configurations preserve the ground state energy of (3n − 4) while making the rewired graphs impossible to optimise with the HT algorithm in the lowgain limit, even for the smallest problem sizes.For the Möbius ladder with no frustration (odd n), the edges J 12 , J N−2,N−3 could be rewired as J 1,N−3 , J 2,N−2 to violate the OSC for any N ≥ 10.Although satisfying the OSC is sufficient for the certain graph structure to be simple, its violation does not necessarily make the instance hard to solve.Other optimisation approaches have to be tested to estimate the relative hardness.
One way to address the relative complexity is to use exact solvers.For example, the commercial solver Gurobi 48 employs various pre-processing techniques and uses heuristics for accelerating the branch-and-bound algorithm 49 that can be applied to mixed-integer programming problems.For problems that cannot be exactly solved for a given time limit, Gurobi evaluates the optimality gap defined as where E best and E lower bound are the best objective and the lower objective bound, respectively.The size of the optimality gap or the time to reach a particular gap could be used as a performance metric for the problem complexity 50 .Hence, the relative hardness of the Ising model on the rewired Möbius ladder graphs can be evaluated by the time it takes Gurobi to reach a zero optimality gap.Another way to evaluate the relative computational hardness is to use heuristic solvers.We consider two physics-inspired algorithms, namely parallel tempering (PT) and chaotic amplitude control.The former [51][52][53][54] is arguably the current state-of-theart algorithm that generally shows better or similar performance over other heuristic methods [55][56][57] , while the latter is inspired by the operation of the Ising machine based on optical parametric oscillators 44,58 .The operational principles of both algorithms and their optimal parameters are discussed in the Methods section.The standard quantitative measure of performance of such stochastic algorithms is the time to solution metric 57,59 , which reflects the time it takes to find the ground state with 99% confidence: where t run is the time for running an algorithm once, R 99 is the number of runs for finding the ground state energy with a probability of 99%, and p gs is the ground state success probability.The computational complexity of the Ising model on 3-regular graphs is interpolated between P and NP classes in Fig. 2c, d, where we use the branch-and-bound, chaotic amplitude control, and PT methods.For the branch-and-bound algorithm within Gurobi solver, Fig. 2c shows the time to zero optimality gap dependence on the percentage of rewired edges in the Möbius ladder graphs of size N ∈ {100, 200, 300}.For all sizes, the initial exponential increase in time is followed by a plateau starting at about 40-50% of rewired edges.For this percentage of rearranged edges, the still recognisable original four-band structure of the Möbius ladder graph has equivalent complexity of random 3-regular graphs for Ising model minimisation.Such equivalence can be associated with frustrated (unsatisfied) edges, namely edges with different signs of s i s j and J ij , the number of which is necessarily minimised at the ground state.Rewiring 40% edges in the Möbius ladder for N = 100 introduces about 8% of frustrated edges, which makes its complexity relatively similar to random 3-regular graphs with 8.6% of frustrated edges.For heuristic algorithms, the optimised time to solution metric is shown in Fig. 2d.For each percentage of rewired edges, we minimise the time to solution with the machine-learner online optimisation package 60 .Unlike the branch-and-bound exact solver, the smaller number of rewired edges is required for achieving a time plateau for the considered stochastic algorithms: the relative hardness of the Möbius ladders with around 2% of rewired edges is equivalent to the complexity of random 3-regular graphs.Similar to Gurobi solver, the computational effort scales exponentially with the number of rewired edges for the chaotic amplitude control method, although a more modest speedup of 2-2.5 times is observed with the PT.
Evidently, all three considered optimisation techniques can successfully determine the relative easiness of the Möbius ladders for the Ising model.The number of rearranged edges in the Möbius ladder graphs, required for achieving the equivalent computational complexity of random 3-regular MaxCut, depends on the method's operational principles.As assessed with arguably state-of-the-art heuristics, this relative complexity of rewired graphs leads to a practically significant result.One may consider an existing physical platform that was previously tested on Möbius ladders and construct the graphs with 2-5% of rewired edges with minimal engineering adjustments.The global Ising minimisation of such rewired Möbius ladder graphs would suggest the physical platform's ability to go beyond the eigenvalue maximisation principle and potentially solve problems that are as hard as random 3-regular MaxCut.Without rewiring, the average computational complexity of the Ising model grows polynomially with the Möbius ladder problem size, as demonstrated for unweighted and weighted graphs with coupling values taken from bimodal and Gaussian distributions in Fig. 2e.

General applicability of the optimisation simplicity criterion.
Any instance of a problem from the P-class is polynomially easy to optimise, while for an arbitrary instance of NP-hard problem, there is no guarantee that the instance is hard.Hardness cannot be guaranteed by violating the proposed OSC, which in itself can only help detect naturally easy instances of NP-hard problems.With the addition of the rewiring procedure proposed above, the relative computational complexity of random graphs can be probed.Till now, the identified simple instances of Ising models were limited to the Möbius ladder graphs.To emphasise the general applicability of the OSC to instances of any NP-hard problem, we show examples of simple graphs in a diverse set of problems that are often chosen to evaluate the performance of the Ising physical machines and computational algorithms.
We apply the OSC to the Ising models with dense, e.g., the Sherrington-Kirkpatrick and the Mattis models, and sparse coupling matrices, where besides 3-regular MaxCut, we examine spin glass models of various topologies including torus, Chimera graph, and 3-regular planar graphs.Where appropriate, in addition to unweighted coupling matrices, we consider commonly chosen probability distributions for interaction strengths such as bimodal, when couplings take values from {−1, 1} with equal probability, and Gaussian when couplings are distributed around zero mean with unit variance (for the full model descriptions, please see Methods).Some of these models belong to the P-class with all instances satisfying the OSC, e.g., the Mattis spin model, unweighted spin glass on a torus, unweighted biased ferromagnet on the Chimera graph, or unweighted ladder graphs with a magnetic field (see Fig. 3).For other models, there exist high chances of getting easy to optimise small size random instances.Across all models, consistently greater probabilities of simple Ising instances are observed for the coupling matrices with values from bimodal and unweighted distributions compared to the Gaussian distribution.Note that in the case of the weighted Möbius ladder graphs, the instances not satisfying the OSC remain easy in terms of both mathematical and computational complexities, as we previously argued.
When testing small-scale Ising simulators, the existence of large fractions of easy instances of NP-hard problems should be taken into account to avoid a misleading assessment of the optimisation capabilities of the platform.A hard random instance would possibly be generated for large problem sizes, while smallscale simulators would likely face easy instances satisfying the OSC.As Fig. 3 shows, the percentage of frustrated edges in the ground state covers the entire range of possible values confirming that the OSC could help identify simple graphs in low and highly frustrated models.

Discussion
Classical and quantum physical systems as analogue simulators can become a superior computational paradigm for solving challenging optimisation problems.Identifying nontrivial and truly hard instances of problems can help evaluate and generalise small size Ising machines' performance.Generally, whether a problem could be considered easy to optimise depends on the existence of an insight into its inherent structure.If there is a way to move through the exponentially large space of possible solutions to the global minimum in polynomial time, then the problem is in the P complexity class.To identify computationally simple instances within the Ising model, we present an optimisation simplicity criterion that is compact and simple to try: one needs to confirm that the signs of the eigenvector, corresponding to the largest eigenvalue of the coupling matrix, coincide with the ground state spin configuration of the Ising model.For instances satisfying the proposed criterion, there is an efficient polynomial time algorithm, e.g., the HT algorithm.For Ising machines employing this algorithm, the quadratic increase in the number of iterations for optimising the Ising model on the Möbius ladder graphs could be expected.
The identification of to optimise Ising coupling matrices allows one to study the continuous computational complexity transition within the same kind of NP-hard problem.For example, the complexity was interpolated between P and NPhard classes for the K-satisfiability problem 61 .To probe an intermediate computational complexity of the maximum cut problem, we introduce a rewiring procedure from the Möbius ladder graphs, the unweighted Ising model on which is in P, to random 3-regular graphs, on which the Ising model is NP-hard.In the case of exact optimisers such as the Gurobi solver, the relative complexity can be evaluated by the size of the optimality gap and time for reaching it, while for heuristic solvers, the time to solution metric could be estimated.For the physical Ising machines, the complexity could be expected to increase exponentially until a certain percentage of rewired edges is reached, making the relative hardness of the rewired graphs similar to random 3-regular graphs.The threshold of the number of rewired edges would depend on a particular machine and its operational principles.
By further exploiting the Gurobi solver, we confirm the superlinear scaling for the average computational complexity of the Ising model on the Möbius ladder graphs with sizes up to 500,000 with couplings from the Gaussian, bimodal and unweighted distributions.The mathematical complexity is also shown to be in P-class for the weighted Ising model on Möbius ladder graphs.
An analysis of the average-case hardness is critical to establishing optimisation performance of small size unconventional computing systems.The general applicability of the proposed criterion for detecting easy to optimise instances of NP-hard problems is demonstrated on a diverse set of the Ising models that includes sparse and dense interaction matrices, weighted and unweighted models, bimodal and Gaussian coupling distributions, with and without a magnetic field, planar and nonplanar geometrical topologies, low and highly frustrated models, regular and not regular graphs.The established fractions of simple instances for such computationally hard problems as 3-regular MaxCut and various spin glass models can be viewed as a lower bound of expected superior optimisation performance for unconventional computing machines on such graphs.The reported simplicity criterion is sufficient but not necessary for an instance to be counted as easy to optimise.Hence, there exist great opportunities for developing other simplicity criteria for identifying easy instances of NP-hard problems.We anticipate that our work will stimulate further studies of the average instance hardness of NP-hard problems and will be followed by other simplicity criteria.
The evidence we provided for the relative computational hardness of specific Ising coupling matrices points to a promising direction for many platforms to reveal their optimisation capabilities to solve complex combinatorial problems.Performance on easy instances of NP-hard problems, satisfying the proposed OSC, does not demonstrate the overall potential of the platform to optimise hard problems and could only confirm the ability of a system to follow one of the largest eigenvectors.This applies to both analogue classical Ising machines and quantum computers, on which quantum-enhanced Ising optimisation has been recently reported 24 .
Selection of the nontrivial or, even better, the hardest instances available in the NP-hard complexity class could tell more about the general optimisation capabilities of physical machines, even of small size, and could lead to more accurate prediction of their large scale performance.Using the rewired Möbius ladder graphs, the optical, electronic and photonic systems could demonstrate their ability to operate beyond the Hopfield-like networks in the low-gain limit and uncover the beneficial internal physical Fraction of instances, satisfying the optimisation simplicity criterion, is shown as a function of problem size N for Gaussian, bimodal, and unweighted coupling distributions.The considered Ising models include Sherrington-Kirkpatrick, 3-regular maximum cut, Mattis spin glass, spin glass on a torus, Möbius ladder graphs, biased ferromagnet on Chimera graph, planar spin glass within a magnetic field.The red dashed line represents models which are polynomially easy to optimise across all problem sizes.For each model, 1000 random matrices are generated per each size, and the ground states are verified with the exact Gurobi solver.
processes for solving challenging problems.As a result, architectures with better optimisation potential will mature faster, approaching the demonstration of computational supremacy on combinatorial optimisation tasks.

Methods
Algorithms. 1. Hopfield-tank neural networks.The numerical integration of the HT algorithm (2) is performed in Fig. 1 by the Euler scheme with the discrete-time step dt = 0.9 and random initial conditions.In all numerical simulations, a hyperbolic tangent is used as an activation function g ¼ tanh x=x 0 À Á and the numerical parameters are τ = 1, I b = 0, x 0 = 3.The polynomial fits are 0.006x 1.986 , 0.01x 1.993 , 0.026x 2.006 for ground state probabilities 50%, 75% and 100%, respectively.We note that for coupling matrices not satisfying the OSC, the HT neural networks could still find energies lower than the energy corresponding to the largest eigenvector by violating the low-gain limit.
2. Gurobi solver.The optimality gaps and times to reach them are obtained with the Gurobi solver on a single core of Intel(R) Core(TM) i9-8950HK CPU 2.0GHz in Fig. 2.
3. Chaotic amplitude control method.The chaotic amplitude control method is an iterative scheme that is inspired by the operation of optical parametric oscillators 44 and based on time evolution of equations: where the signs of x i represent the Ising spins, p is the linear gain, ϵ is the coupling strength scaling coefficient, e i is the error signal, β is the rate of change of error signal, and a is the target amplitude.The dynamics of x i elements in Eq. ( 9) is similar to a Hopfield network evolution with the addition of nonlinearity, while the amplitude alignment of the network elements as in Eq. ( 10) was first argued to be necessary for minimisation of both discrete and continuous spin Hamiltonians for gain-dissipative simulators in 62 .Following 44 , the parameters a, p, and β are dynamically adjusted as where the baseline of the target amplitude a is set to one, δ is the sensitivity to energy variations, ΔE = E best − E(t) is the difference between the best-found Ising energy and the current energy at time t, p 0 is the linear gain baseline.The parameter β increases with a positive rate γ for the maximum allowed time τ, otherwise is set to zero and t c is set to t, where t c is the last time when the bestknown energy E best was updated or β was reset.The numerical simulation of Eqs. ( 9)-( 13) could efficiently sample the low energy states of Ising Hamiltonian 44 and was recently implemented on a field-programmable gate array as chaotic amplitude control method 58 .We note that the parameters e i in Eq. ( 9) play a critical role in the algorithm's performance and make the local minima escape mechanism of the chaotic amplitude control method somewhat similar to that of PT.Namely, the error signals e i could be seen as a set of inverse temperatures, each applied to individual spin.Hence, the energy landscape exploration within subspaces is performed during network elements' evolution at different temperatures.For spins whose amplitude x 2 i is much less than the target amplitude a, the signals e i are large and help the system to quickly settle in a local energy minimum within the phase space of these spins (small temperature regime with a rough energy landscape for a subset of spins).For spins with amplitudes close to the target value, the signals e i are small and facilitate crossings over energy barriers within this phase subspace (high temperature regime with a smooth energy landscape for a subset of spins).The global optimisation could be then achieved due to fluctuations of 'inverse temperatures' e i in time, which realise a continuous exploration of random spin subspaces at low and high temperatures.The importance of multiplying the interactions by amplitudes e i has also been argued in terms of the destabilisation of low energy local minima 58 and due to the energy-conserving rotary motion 63 .
4. Parallel tempering.PT, or replica exchange Monte Carlo method, takes advantage of the dependence of the energy landscape on the temperature of the system: a rough landscape at low temperatures becomes smooth at high enough temperatures.In this algorithm [51][52][53][54] , multiple states (replicas) of a particular Ising model are simulated independently at different temperatures with a Monte Carlo algorithm.At low temperatures, the local exploration of an energy landscape is realised, and replicas could quickly get trapped in local minima.In contrast, a global exploration of energy profile is facilitated with replicas easily crossing energy barriers at high temperatures.The global optimisation could be then achieved by exchanging replicas at low and high temperatures.For the optimal performance of PT, one needs to ensure the exchange in the reverse direction, which requires a trade-off between two factors.On the one hand, the exchange mechanism is accelerated once the acceptance probabilities for exchanges are high, which can usually be achieved by increasing the number of replicas.On the other hand, the larger number of replicas slows down the algorithm and increases the exchange time between replicas at the lowest and highest temperatures.
Optimal parameters.For PT and chaotic amplitude control methods, the optimal sets of parameters are determined with the machine-learner online optimisation package (m-loop) 60 over ten random problem instances for each rewiring percentage per each problem size which results in a total of more than 150 graphs optimised for each problem size.The time to solution has been used as a target function for optimisation with the success probability calculated by running algorithms 30 times for each instance.The number of m-loop steps was fixed to 200.Both algorithms are implemented in Python and translated to optimised machine code with Numba, while Numba-compiled algorithms are known to approach a performance similar to Fortran.The reported time to solution performance is achieved on a single core of Intel(R) Core(TM) i9-8950HK CPU 2.0GHz.
For chaotic amplitude control method, the optimal time to solution for a given rewiring percentage for each problem size is determined by optimising the following seven parameters with the m-loop: number of iterations (N iter ), maximum number of iterations without energy change (N τ ), coupling strength (ϵ), linear gain baseline (p 0 ), rate of increase of β (γ), sensitivity to energy variations (δ), and time step (dt).The averaged optimised parameters are listed in Table 1, where the linear gain baseline is determined through the parameter α and the largest eigenvalue λ 0 of coupling matrix J as p 0 = 1 − αλ 0 .
For PT, the optimal time to solution for a given rewiring percentage for each problem size is determined by optimising the following five parameters with the mloop: the number of Monte Carlo sweeps (N MC ), number of replicas (N R ), number of iterations for equilibration (N eq ), low temperature (T low ) and high temperature (T high ).Temperatures at each replica are set based on the geometric schedule 54,57 .Compared to the pseudocode in 57 , our PT implementation includes an additional extra hyperparameter N eq , which allows individual replicas to equilibrate over a few Monte Carlo sweeps before exchanging states of neighbouring pairs of replicas.The Metropolis update mechanism is used for the Monte Carlo simulations.The average values of optimal parameters are shown in Table 2.
Ising models.In Fig. 3, the non-exhaustive list of problems in which one can find polynomially easy Ising instances includes: 1. Sherrington-Kirkpatrick (SK) model of spin glasses 64 .The fully-connected SK instances have a coupling matrix with elements from Gaussian distribution with zero mean and unit variance (Gaussian-SK).The Gaussian-SK model is NP-hard 65 The individual sets of parameters are selected for each rewiring of 3-regular graphs in Fig. 2 by optimising time to solution for ten rewired graphs with the m-loop method.The presented values are averaged over the percentage of rewired edges for each problem size.The optimised parameters include number of iterations (N iter ), the maximum number of iterations without energy change (N τ ), coupling strength (ϵ), coefficient for the linear gain baseline (α), rate of increase of change of error signal (γ), sensitivity to energy variations (δ), and time step (dt).The individual sets of parameters are selected for each rewiring of 3-regular graphs in Fig. 2 by optimising time to solution for ten rewired graphs with the m-loop method.The presented values are averaged over the percentage of rewired edges for each problem size.The optimised parameters include the number of Monte Carlo sweeps (N MC ), number of replicas (N R ), number of iterations for equilibration (N eq ), low temperature (T low ) and high temperature (T high ).
though the ground state with precision of − δ) can be found in polynomial time for any δ > 0 when the coupling coefficients are taken from the Gaussian distribution with zero mean and variance σ = 1/N 66 .The probability of finding an easy instance of the Gaussian-SK problem with the OSC decreases from 45-100% for size N = 3-10 to 10-20% for 20-25 size.The SK model stays in the NP-Hard class 67 when the coupling values are chosen from bimodal distribution (bimodal-SK).In this case, the probability of easy instances drops from 65-100% to 20% for problem sizes 3-10 and 20-25, respectively.Both models have 100% simple instances for N = 3 and all instances are simple for N = 5 in the case of bimodal distribution.The unweighted SK model coincides with the complete unweighted graphs, which were considered for the complexity continuum transition of kregular graphs and argued to be polynomially simple.We note that the ground states of complete graphs of odd size starting from N = 43 can be confirmed up to 1 frustrated edge with Gurobi solver in 1200 s, so they were additionally verified with the recent physics-inspired algorithms 44,62 .Both Gaussian-SK and bimodal-SK are commonly chosen for comparing physical machines 22,68 and computational algorithms 44,57,69 .
2. Mattis spin glass (Mattis SG) model 70 .In the Mattis model, random variables ϵ are generated for each site i according to a specified probability distribution to build separable spin interactions as J ij = f(R ij )ϵ i ϵ j , where f(R ij ) is the adjacency matrix that specifies the topology of a graph.Such a model does not have frustrations, and the ground state is identical to the configuration of the random variables s i = ϵ i .Also, one may notice that the Mattis model is equivalent to gauge transformation J gauged ij ¼ J F ij ϵ i ϵ j which conceals the planted ground state of the problem with ferromagnetic couplings J F ij .For both Gaussian and bimodal probability distributions of couplings, all Mattis spin models' instances satisfy the OSC, which generalises to any problem size.Thus, the Mattis SG belongs to the Pclass.The Mattis model was recently used for evaluating the performance of photonic Ising machines 11,68 .
3. Maximum cut on 3-regular graphs.In addition to unweighted 3-regular graphs, we considered 3-regular MaxCut with couplings from bimodal and Gaussian distributions.The bimodal 3-regular MaxCut exhibits a similar probability of easy instances as unweighted 3-regular graphs, while the probabilities for Gaussian 3-regular MaxCut are slightly higher on average than for Gaussian-SK.Besides, the case of 3-regular graphs on the Möbius ladder is considered for bimodal and Gaussian coupling distributions.The MaxCut problems are commonly chosen for evaluating physical simulators 21,22,24,30,71 .
4. Spin glass model on a torus (SG-torus).A torus is represented by a twodimensional rectangular lattice with periodic boundaries in both directions and nearest-neighbour interactions.The unweighted SG-torus model satisfies the OSC for any problem size.The Gaussian SG-torus is less likely to have simple graphs compared to Gaussian-SK, while the chances of about 40% hold even for a problem size of N = 40 for bimodal SG-torus.The SG-torus models were recently used for comparing the large scale performance of optimisation physics-inspired algorithms 57 .
5. Planar spin glass within a magnetic field.One of the earliest proofs of NPhardness of the Ising model was demonstrated for a three-dimensional spin glass and a planar spin glass within a uniform magnetic field h i = −1 and unweighted antiferromagnetic interactions 18 .Conveniently for us, the Möbius ladder graphs can be easily rewired to planar cubic graphs by avoiding the twist and becoming ladder graphs.All unweighted ladder graphs with a magnetic field satisfy the OSC.By exploiting the rewiring procedure with an additional planarity constraint, about 50% of random planar 3-regular graphs happen to be simple for a problem size of 20.We also note that all found planar graphs of size six are simple graphs.
6. Biased ferromagnet on Chimera graph (BF-Chimera).The model represents an unweighted ferromagnetic coupling matrix on Chimera graph with fields p(h i = 0) = p 0 and p(h i = 1) = p 1 where p 0 ≫ p 1 that bias s i = 1 for all spins as the global optimal solution.This model was introduced as a toy example to get an intuition behind the optimisation behaviour of the D-Wave machine and classical algorithms 50 .The BF-Chimera model has no frustration and its instances satisfy the OSC and thus are in P-class.Though this is the only model in our list that was not argued to be hard before, its presence here could serve for studying the complexity of other known Ising models with Chimera topology.We note that there could be an additional overhead due to topological embedding in simulators where the long-range interactions are nontrivial to engineer, e.g., the Chimera graphs for some Ising problems.

Fig. 1
Fig. 1 Optimisation of the unweighted Ising model on the Möbius ladder graphs with the Hopfield-Tank algorithm.a Illustration of the Möbius ladder graph on Möbius strip and on the circular graph.Two possible frustrated edges in the ground state are highlighted in red.b The evolution of amplitudes v i and the corresponding proximity to the ground state (GS) are shown for the Möbius ladder graph of size N = 1000 over N iter = 3000 time iterations of the Hopfield-Tank algorithm.All low energy levels E λ i correspond to the projected eigenvectors sign(e i ) of the distinct largest eigenvalues λ i .c The number of time iterations N iter of the Hopfield-Tank algorithm is shown for optimising the Möbius ladder graphs of sizes up to N = 10,000 with desired ground state probability ranges of p gs ∈ {50-55%, 75-80%, 99-100%}.The solid lines correspond to a quadratic fit confirming that the Ising model on Möbius ladder graphs can be solved in polynomial time.The number of algorithm runs per each graph size is fixed to 250.d The ground state probabilities as a function of Möbius Ladder size are shown for the fixed number of time iterations N iter ∈ {10,000, 50,000, 250,000}.

Fig. 2
Fig. 2 Continuum complexity transition of the Ising model on 3-regular graphs.a, b The rewiring procedure of two edges for violating the optimisation simplicity criterion in the Möbius ladder graphs of size N = 2n for any even N in (a) and most even N in (b).The removed and added edges are shown with red solid and dashed green lines, respectively.c The Möbius ladder graphs are depicted in the upper panel with 10%, 60% and 100% of rearranged edges.The relative computational hardness of the unweighted Ising model on the rewired Möbius ladder graphs is evaluated by the median time required for reaching zero optimality gap with Gurobi solver and is shown in the lower panel for problem sizes 100, 200 and 300.The 100 random graphs are optimised for each problem size for every percentage of rewired edges with shaded regions indicating an interquartile range.d The median time to solution metric as a function of the percentage of rewired edges in the Möbius ladder graphs is shown for chaotic amplitude control and parallel tempering methods.e The average computational hardness is evaluated by the median time required for reaching zero optimality gap with Gurobi solver for the Möbius ladder graphs of sizes up to 500,000 with couplings from the unweighted, bimodal, and Gaussian distributions.The ten random graphs are optimised for each problem size.

Fig. 3
Fig.3Probability of finding easy instances for various Ising models.Fraction of instances, satisfying the optimisation simplicity criterion, is shown as a function of problem size N for Gaussian, bimodal, and unweighted coupling distributions.The considered Ising models include Sherrington-Kirkpatrick, 3-regular maximum cut, Mattis spin glass, spin glass on a torus, Möbius ladder graphs, biased ferromagnet on Chimera graph, planar spin glass within a magnetic field.The red dashed line represents models which are polynomially easy to optimise across all problem sizes.For each model, 1000 random matrices are generated per each size, and the ground states are verified with the exact Gurobi solver.

Table 1
Optimal parameters for chaotic amplitude control method.

Table 2
Optimal parameters for parallel tempering.