Scaling advantage of nonrelaxational dynamics for high-performance combinatorial optimization

Many complex systems such as spin glasses, interacting proteins, large scale hardware, and financial portfolios, can be described as ensembles of disordered elements that have competing frustrated interactions[1] and rugged energy landscapes. There has been a growing interest in using physical simulators called “Ising machines” in order to reduce time and resources needed to identify configurations that minimize their total interaction energy, notably that of the Ising Hamiltonians $\mathcal{H}$ with $\mathcal{H}(\boldsymbol{\sigma})=-\frac{1}{2}\sum_{ij}\omega_{ij}\sigma_{i}\sigma_{j}$ (with $\omega_{ij}$ the symmetric Ising couplings, i.e., $\omega_{ij}=\omega_{ji}$, and $\sigma_{i}=\pm 1$) that is related to many nondeterministic polynomial-time hard (NP-hard) combinatorial optimization problems and various real-world applications[2] (see [3] for a review). Recently proposed implementations include memresistor networks[4], micro- or nano-electromechanical systems[5], micro-magnets[6, 7], coherent optical systems[8], hybrid opto-electronic hardware[5, 10, 11], integrated photonics[12, 13, 14], flux qubits[15], and Bose-Einstein condensates[16]. In principle, these physical systems often possess unique properties, such as coherent superposition in flux qubits[17] and energy efficiency of memresistors[18, 4], which could lead to a distinctive advantage compared to conventional computers (see Fig. 1 (a)) for the sampling of low energy states. In practice, the difficulty in constructing connections between constituting elements of the hardware is often the main limiting factor to scalability and performance for these systems[19, 15]. Moreover, these devices often implement schemes that are directly related to the concept of annealing (either simulated[20, 21], mean-field[2, 3], chaotic[18, 24], and quantum[25, 17]) in which the escape from the numerous local minima[26] and saddle points[27] of the free energy function can only be achieved under very slow modulation of the control parameter (see Fig. 1 (b)). These methods are dependent on non-equilibrium dynamics called aging that, according to recent numerical studies[28], is strongly non-ergodic and seems to explore only a confined subspace determined by the initial condition rather than wander in the entire configurational space[29] for mean-field spin glass models. In general, such systems find the solutions of minimal energy only after many repetitions of the relaxation process.

Interestingly, alternative dynamics that is not based on the concepts of annealing and relaxation may perform better for solving hard combinatorial optimization problems[30, 31, 32]. Various kinds of dynamics have been proposed[33, 34, 35, 36, 3], notably chaotic dynamics[37, 38, 18, 39, 40], but have either not been implemented onto specialized hardware[41, 37] or use chaotic dynamics merely as a replacement to random fluctuations[38, 18]. We have recently proposed that the control of amplitude in mean-field dynamics can significantly improve the performance of Ising machines by introducing error correction terms (see Fig. 1 (d)), effectively doubling the dimensionality of the system, whose role is to correct the amplitude heterogeneity[30]. Because of the similarity of such dynamics with that of a neural network, it can be implemented especially efficiently in electronic neuromorphic hardware where memory is distributed with the processing[42, 43, 44]. In this paper, we show that this nonrelaxational dynamics is able to escape at a much faster rate than relaxational ones from local minima and saddles of the corresponding energy function. Importantly, the exponential scaling factor with respect to system size of the time needed to reach ground-state configurations of spin glasses, called time to solution, is significantly smaller in the former case, which raises the question whether this nonrelaxational dynamics is qualitatively different from the very slow relaxation observed in classic Monte Carlo simulations of spin glasses. In order to extend numerical analysis to large problem sizes and limit finite-size effects, we implement a version of nonrelaxational dynamics that we name chaotic amplitude control (CAC) on a field programmable gate array (FPGA, see Fig. 1 (c)) and show that the developed hardware can be faster for finding ground-states in the limit of large problem sizes than many state-of-the-art algorithms and Ising machines for some reference benchmarks with enhanced energy efficiency.

For the sake of simplicity, we consider the classical limit of Ising machines for which the state space is often naturally described by analog variables (i.e., real numbers) noted $x_{i}$ in the following. The variables $x_{i}$ represent measured physical quantities such as voltage[4] or optical field amplitude[8, 5, 10, 11, 12, 13] and these systems can often be simplified to networks of interacting nonlinear elements whose time evolution can be written as follows:

It is well known that, when $\beta_{i}=\beta$ for all $i$ and the noise is not taken into account ($\sigma_{0}=0$), the time evolution of this system is motion in the state space that seeks out minima of a potential function[1] (or Lyapunov function) $V$ given as $V=\beta\mathcal{H}(\boldsymbol{y})+\sum_{i}V_{b}(y_{i})$ where $V_{b}$ is a bistable potential with $V_{b}(y_{i})=-\int_{0}^{y_{i}}f_{i}(g_{i}^{-1}(y))\mathop{}\!\mathrm{d}y$ and $\mathcal{H}(\boldsymbol{y})=-\frac{1}{2}\sum_{ij}\omega_{ij}y_{i}y_{j}$ is the extension of the Ising Hamiltonian in the real space with $y_{i}=g_{i}(x_{i})$ (see Supplementary Materials S1.1). The bifurcation parameter $\beta$, which can be interpreted as the inverse temperature of the naive TAP equations[48], the steepness of the neuronal transfer function in Hopfield-Tank neural networks[45], or to the coupling strength in coherent Ising machines[8, 5, 10], is usually decreased gradually in order to improve the quality of solutions found. This procedure has been called mean-field annealing[3], and can be interpreted as a quasi-static deformation of the potential function $V$ (see Fig. 1 (b)). There is, however, no guarantee that a sufficiently slow deformation of the landscape $V$ will ensure convergence to the lowest energy state contrarily to the quantum adiabatic theorem[50] or the convergence theorem of simulated annealing[51]²²2At fixed $\beta$, global convergence to the minimum of the potential $V$ can be assured if $\sigma_{0}$ is gradually decrease with $\sigma_{0}(t)^{2}\sim\frac{c}{log(2+t)}$ and $c$ sufficiently large[52]. The parameter $\sigma_{0}^{2}$ is analogous to the temperature in simulated annealing in this case. The global minimum of the potential $V$ does not, however, generally correspond to that of the Ising Hamiltonians $\mathcal{H}$ at finite $\beta$.. Moreover, the statistical analysis of spin glasses suggests that the potential $V$ is highly non-convex at low temperature and that simple gradient descent very unlikely reaches the global minimum of $\mathcal{H}(\boldsymbol{\sigma})$ because of the presence of exponentially numerous local minima[26] and saddle points[27] as the size of the system increases. The slow relaxation time of Monte-Carlo simulations of spin glasses, such as when using simulated annealing, might also be explained by similar trapping dynamics during the descent of the free energy landscape obtained from the TAP equations[27]. In the following, we consider in particular the soft spin description obtained by taking $f_{i}(x_{i})=(-1+p)x_{i}-x_{i}^{3}$ and $y_{i}=g(x_{i})=x_{i}$, where $p$ is the gain parameter, which is the canonical model of the system described in eq. (1) at proximity of a pitchfork bifurcation with respect to the parameter $p$. In this case, the potential function $V_{b}$ is given as $V_{b}(x_{i})=(-1+p)\frac{x_{i}^{2}}{2}-\frac{x_{i}^{4}}{4}$ and eq. (1) can be written as $\frac{\mathop{}\!\mathrm{d}x_{i}}{\mathop{}\!\mathrm{d}t}=-\frac{\partial V}{\partial x_{i}}$, $\forall i$.

In order to define nonrelaxational dynamics that is inclined to visit spin configurations associated with lower Ising Hamiltonian, we introduce error signals, noted $e_{i}\in\mathbb{R}$, that modulate the strength of coupling $\beta_{i}$ to the $i$th nonlinear element such that $\beta_{i}(t)$ defined in eq. (1) is expressed as $\beta_{i}(t)=\beta e_{i}(t)$ with $\beta>0$. The time evolution of the error signals $e_{i}$ are given as follows[30]:

The dynamics described by the coupled equations (1) and (2) is not derived from a potential function because error signals $e_{i}$ introduce asymmetric interactions between the $x_{i}$ and the computational principle is not related to a gradient descent. Rather, the addition of error variables results in additional dimensions in the phase space via which the dynamics can escape local minima. The mechanism of this escape can be summarized as follows. It can be shown (see the Supplementary Materials S1.2) that the dimension of the unstable manifold at equilibrium points corresponding to local minima $\boldsymbol{\sigma}$ of the Ising Hamiltonian depends on the number of eigenvalues $\mu(\boldsymbol{\sigma})$ with $\mu(\boldsymbol{\sigma})>F(a)$ where $\mu(\boldsymbol{\sigma})$ are the eigenvalues of the matrix $\{\frac{\omega_{ij}}{|h_{i}|}\}_{ij}$ (with internal field $h_{i}$) and $F$ a function given as $F(y)=\frac{\psi^{\prime}(y)}{\psi(y)}y$ and $\psi(y)=\frac{f(g^{-1}(y))}{(g^{-1})^{\prime}(y)}$. Thus, there exists a value of $a$ such that all local minima (including the ground state) are unstable and for which the system exhibits chaotic dynamics that explores successively candidate boolean configurations. The energy is evaluated at each step and the best configuration visited is kept as the solution of a run. Interestingly, this chaotic search is particularly efficient for sampling configurations of the Ising Hamiltonian close to that of the ground state using a single run although the distribution of sampled states is not necessarily given by the Boltzmann distribution. Note that the use of chaotic dynamics for solving Ising problems has been discussed previously[53, 18], notably in the context of neural networks, and it has been argued that chaotic fluctuations may possess better properties than Brownian noise for escaping from local minima traps. In the case of the proposed scheme, the chaotic dynamics is not merely used as a replacement to noise. Rather, the interaction between nonlinear gain and error-correction results in the destabilization of states associated with lower Ising Hamiltonian.

Ensuring that fixed points are locally unstable does not guarantee that the system does not relax to periodic and chaotic attractors. We have previously proposed that non-trivial attractors can also be destabilized by ensuring the positive rate of entropy production using a modulation of the target amplitude[30]. In this paper, we propose an alternative heuristic modulation of the target amplitude that is better suited for a digital implementation than the one proposed in [30]. Because the value of $a$ for which all local minima is unstable is not known a priori, we propose instead to destabilize the local minima traps by dynamically modulating $a$ depending on the visited configurations $\boldsymbol{\sigma}$ as follows:

In order to verify that the nonrelaxational dynamics of chaotic amplitude control is able to accelerate the search of mean-field dynamics for finding the ground state of typical frustrated systems, we look for the ground-states of Sherrington-Kirkpatrick (SK) spin glass instances using the numerical simulation of eqs. (1) to (3) and compare time to solutions with those of two closely related relaxational schemes: noisy mean-field annealing (NMFA)[2] and the simulation of the coherent Ising machine (simCIM) ³³3See Supplementary Materials S4.1-2. Because the arithmetic complexity of calculating one step of these three schemes is dominated by the matrix-vector multiplication (MVM), it is sufficient for the sake of comparison to count the number of MVM, noted $\nu$, to find the ground state energy of a given instance with $99\%$ success probability, with $\nu(K)=K\frac{\text{ln}(1-0.99)}{\text{ln}(1-p_{0}(K))}$ and $p_{0}(K)$ the probability of visiting a ground state configuration at least once after a number $K$ of MVMs in a single run. In Fig. 2, NMFA (a) and the CAC (b) are compared using the averaged success probability $\langle p_{0}\rangle$ of finding the ground state for 100 randomly generated SK spin glass instances per problem size $N$. Note that the success probability of the mean-field annealing method does not seem to converge to 1 even for large annealing time (see Fig. 2 (a)). Because the success probability of NMFA and simCIM remains low at larger problem size, its correct estimation requires simulating a larger number of runs which we achieved by using GPU implementations of these methods. On the other hand, the average success probability $\langle p_{0}\rangle$ of CAC is of order 1 when the maximal number of MVM is large enough, suggesting that the system rarely gets trapped in local minima of the Ising Hamiltonian or non-trivial attractors. In Figure 2 (c) and (d) are shown the $q^{\rm th}$ percentile (with $q=50$, i.e., the median) of the MVM to solution distribution, noted $\nu_{q}(K;N)$, for various duration of simulation $K$, where $K$ is the number of MVMs of a single run. The minimum of these curves, noted $\nu_{q}^{*}(N)$ with $\nu_{q}^{*}(N)=\text{min}_{K}\nu_{q}(K;N)$, represents the optimal scaling of MVM to solution vs. problem size $N$[15]. Using the hypothesis of an exponential scaling with the square root of problem size $N$, CAC exhibits significantly smaller scaling exponent ($\gamma=0.18\pm 0.06$) than the NMFA ($\gamma=0.47\pm 0.04$) and simCIM ($\gamma=0.54\pm 0.03$, see inset in Fig. 2 (e)). We have verified that this scaling advantage holds for various parameters of the mean-field annealing (see Supplementary Materials S4.1). Note that a root-exponential scaling of the median time to solution has been previously reported for SK spin glass problems[15, 54] and other NP-Hard problems[55].

Although comparison of algebraic complexity indicates that CAC has a scaling advantage over mean-field annealing, it is in practice necessary to compare its physical implementation against other state-of-the-art methods because the performance of hardware depends on other factors such as memory access and information propagation delays. To this end, CAC is implemented into a FPGA because its similarity with neural networks makes it well-fitted for a design where memory is distributed with processing (see Supplementary Materials S2 for the details of the FPGA implementation). The organization of the electronic circuit can be understood using the following analogy. Pairs of analog values $x_{i}$ and $e_{i}$, which represent averaged activity of two types of neurons, are encoded within neighboring circuits. This micro-structure is repeated $N$ times on the whole surface of the chip which resembles the columnar organization of the brain. The nonlinear processes $f_{i}(x_{i})$, which model the local-population activation functions and are independent for $i\neq j$, are calculated in parallel. The coupling between elements $i$ and $j\in\{1,\ldots,N\}$ that is achieved by the dot product in eq. (1) is implemented by circuits that are at the periphery of the chip and are organized in a summation tree reminiscent of dendritic branching (see Fig. 1 (c)). The power consumption of the developed hardware never exceeds 5W because of limitations of the development board that we have used.

First, we compare the FPGA implementation of CAC against state-of-the-art CPU algorithms: break-out local search[56] (BLS) that has been used to find many of the best known maximum-cuts (equivalently, Ising energies) from the GSET benchmark set[57], a well-optimized single-core CPU implementation of parallel tempering (or random replica exchange Monte-Carlo Markov chain sampling)[58, 59] (PT, courtesy of S. Mandra), simulated annealing (SA)[60]. Figure 3 (a) shows that the CAC on FPGA has the smallest real time to solution $\tau_{q}^{*}$ against most other state-of-the-art algorithms despite just 5W power-consumption where $\tau_{q}^{*}(N)$ is the optimal $q^{\rm th}$ percentile of time to solution with 99% success probability and is given as $\tau_{q}^{*}(N)=\text{min}_{T}\tau_{q}(T;N)$ where $\tau(T)$ of a given instance is $\tau(T)=T\frac{\text{ln}(1-0.99)}{\text{ln}(1-p_{0}(T))}$ and $T$ is the duration in seconds of a run. The probability $p_{0}(T)$ is evaluated using 100 runs per instance. The results of the CAC algorithm run on a CPU are also included in Fig. 3. The CPU implementation of CAC written in python for this work is not optimized and is consequently slower than other algorithms. However, its scaling of time to solution with problem size is consistent with that of CAC on FPGA. Figure 3 shows that CAC implemented on either CPU or FPGA has a significantly smaller increase of time to solution with problem size than SA run on CPU. Note that the power consumption of transistors in the FPGA and CPU scales proportionally to their clock frequencies. In order to compare different hardware despite the heterogeneity in their power consumption, the $q^{\rm th}$ percentile of energy-to-solution $E_{q}^{*}$, i.e., the energy $E_{q}^{*}$ required to solve SK instances with $E_{q}^{*}=P\tau_{q}^{*}$ and $P$ the power consumption⁴⁴4For the sake of simplicity, we assume a 20 watts power consumptions for the CPU. These numbers represent typical orders to magnitude for contemporary digital systems., is plotted in Fig. 3 (b). CAC on FPGA is $10^{2}$ to $10^{3}$ times more energy efficient than state-of-the-art algorithms running on classical computers.

The Monte Carlo methods SA and PT have moreover been recently implemented on a special-purpose electronic chip called Digital Annealer (DA)[20]. In Fig. 4, we show the scaling exponents of 50th and 80th percentiles of the time to solution distribution for problem sizes $N=800$ to $N=1100$ based on the hypothesis of scaling in $e^{\gamma N}$ obtained by fitting data shown in Figs. 2 and 3 (see “CAC fully parallel” and ‘CAC 100$\times$100”, respectively, in Fig. 4 for the numerical values of the scaling exponents⁵⁵5We replicated the benchmark method of [20] for the sake of the comparison by fitting time to solution from $N=800$ up to $N=1100$.), and compare them to that reported for SA and PT implemented on CPU and DA. The scaling obtained from fitting the time to solution in Fig. 2 is based on the assumption that the matrix-vector multiplication can be calculated fully in parallel in a time that scales as $\text{log}(N)$ instead of $N^{2}$ (see Materials and Methods section) at least up to $N=1100$. We include this hypothesis because many other Ising machines exploit the parallelization of matrix-vector multiplication for speed up[20, 61], whereas the current implementation of CAC iterates on block matrices of size 100 by 100 and is thus only partially parallel because of resource limitations specific to the downscale FPGA used in this work. Note that the time to compute the matrix-vector multiplication is not the dominant term in the exponential scaling behavior of time to solution at large $N$. The scaling of time to solution for the nonrelaxational dynamics observed is significantly smaller than the ones of standard Monte Carlo methods SA and PT[20], especially in the case of a fully parallel implementation. The scaling exponents of fully parallel CAC is smaller than that of DA and on par with that of PT on DA (PTDA), although CAC does not require simulating replica of the system and is thus faster in absolute time than PTDA[20].

Next, the proposed implementation of nonrelaxational dynamics is compared to other recently developed Ising machines (see Fig. 5). The relatively slow increase of time to solution with respect to the number of spins $N$ when solving larger SK problems using CAC suggests that our FPGA implementation is faster than the Hopfield neural network implemented using memresistors (mem-HNN)[4], the restricted Boltzmann machine using a FPGA (FPGA-RBM)[54] at large $N$. Extrapolations are based on the hypotheses of scaling in $e^{\gamma N}$ and $e^{\gamma\sqrt{N}}$ by fitting the available experimental data up to $N=100$ for mem-HNN and FPGA-RBM, $N=150$ for NTT CIM, and $N=1100$ for FPGA-CAC. Figure 5 shows that mem-HNN, FPGA-RBM, and NTT CIM have similar scaling exponents, although FPGA-RBM tends to exhibit a scaling in $e^{\gamma N}$ rather than $e^{\gamma\sqrt{N}}$ for $N\approx 100$[54]. It can be nonetheless expected that the algorithm implemented in mem-HNN, which is similar to mean-field annealing, has the same scaling behavior as simCIM and NMFA (see Fig. 2).

It is noteworthy to mention that a recent implementation of the simulated bifurcation machine[61] (SBM) which is not based on the gradient descent, similarly to CAC, but based on adiabatic evolutions of energy conservative systems performs well in solving SK problems. Both SBM and CAC exhibit smaller time to solution than other gradient based methods. SBM has been implemented on a FPGA (the Intel Stratix 10 GX) that has approximately 5 to 10 times more adaptive logic modules than the KU040 FPGA used to implement CAC. In order to compare SBM and CAC if implemented on an equivalent FPGA hardware, we plot in Fig. 5 the estimation of the time to solution for a fully parallel implementation of CAC using the hypothesis that one matrix-vector multiplication of size $1100\times 1100$ can be achieved in $0.3\mu s$. This is the same time to compute a MVM that we can infer from time to solution reported in [61] for SBM with binary connectivity given that problems of size $N=100$ ($N=700$) are solved in $29\mu s$ ($55ms$) and 94 (81000) MVMs, respectively. Note that SBM can reach the ground-states after approximately 20 times less MVMs than CAC at $N=100$ but only 5 times less MVM at $N=1000$, suggesting that the speed of SBM depends largely on hardware rather than an algorithmic scaling advantage. Moreover, the simulated bifurcation machine[61] does not perform significantly better than our current implementation of CAC for solving instances of the reference MAXCUT benchmark set called GSET[57] (see Tab. 1) even compared to the case of the implementation on the smaller KU040 FPGA with the probability of finding maximum cuts of the GSET in a single run that is much smaller with SBM. Comparison of the scaling behavior of time to solution between CAC and SBM is unfortunately not possible based on available data[61].

Lastly, we consider the whole distribution of time to solution in order to compare the ability of various methods to solve harder instances. As shown in Fig. 6 (a), the cumulative distribution function (CDF) $P(\tau;T)$ of time to solution with $99\%$ success probability $\tau$ is not uniquely defined as it depends on the duration $T$ of the runs. We can define an optimal CDF $P^{*}(\tau)$ that is independent of the runtime $T$ as follows: $P^{*}(\tau)=\text{max}_{T}P(\tau;T)$. Numerical simulations show that this optimal CDF is well described by lognormal distribution, that is $P^{*}(log(\tau))\sim\mathcal{N}(\mu,\sqrt{v})$ where $\sqrt{v}$ is the standard deviation of $log(\tau)$ (see Figs. 6 (b), (c), and (d) for the cases of CAC, SA, and NMFA, respectively). In Fig. 6 (e), it is shown that the scaling of the standard deviation $\sqrt{v}(N)$ with the problem size $N$ is significantly smaller for CAC, which implies that harder instances can be solved relatively more rapidly than using other methods. This result confirms the advantageous scaling of higher percentiles for CAC that was observed in Figs. 2 and 3.

The framework described in this paper can be extended to solve other types of constrained combinatorial optimization problems such as the traveling salesman[45], vehicle routing, and lead optimization problems. Moreover, it can be easily adapted to a variety of recently proposed Ising machines[4, 5, 6, 7, 8, 5, 10, 11, 12, 13, 15] which would benefit from implementing a scheme that does not rely solely on the descent of a potential function. In particular, the performance of CIM[5, 10], mem-HNN[4], and chip-scale photonic Ising machine[14], which have small time to solution for small problem sizes⁶⁶6Relaxation dynamics may be faster for solving small problem sizes for which it may be sufficient to do rapid sampling based on convex optimization[63]. ($N\approx 100$) but with a relatively large scaling exponent that limit their scalability, could be significantly improved by adapting the implementation we propose if these hardware can be shown to be able to simulate larger numbers of spins experimentally. Rapid progress in the growing field of Ising machines may allow to verify scaling behaviors of the various methods at larger problem sizes and, thus, limit further finite-size effects.

The scaling exponents we have reported in this paper are based on the integration of the chaotic amplitude control dynamics using a Euler approximation of its ODEs. We have noted that the scaling of MVMs to solution of SK spin glass problems is reduced when the Euler time step is decreased (see Supplementary Materials S2.8). The scaling exponents of CAC might thus be smaller than reported in this paper in the limit of a more accurate numerical integration over continuous time. It is therefore important future work to evaluate the scaling for $N\gg 1000$ using a faster numerical simulation method. Such numerical calculations require a careful analysis of the integration method of ODEs, numerical precision, and tuning of parameters. It is also of considerable interest to implement CAC on an analog physical system for further reduction of power consumption.

Nonrelaxational dynamics described herein is not limited to artificial simulators and likely also emerge in natural complex systems. In particular, it has been hypothesized that the brain operates out of equilibrium at large scales and produces more entropy when performing physically and cognitively demanding tasks[64, 65] in particular for the ones involving creativity for maximizing reward rather than memory retrieval. Such neural processes cannot be explained simply by the relaxation to fixed point attractors[1], periodic limit cycles[66], or even low-dimensional chaotic attractors whose self-similarity may not be equivalent to complexity but set the conditions for its emergence[67, 65]. Similarly, evolutionary dynamics that is characterized as non-ergodic when the time required for representative genome exploration is longer than available evolutionary time[68] may benefit from nonrelaxational dynamics rather than slow glassy relaxation for faster identification of high-fitness solutions. A detailed analytic comparison between the slow relaxation dynamics observed in Monte Carlo simulations of spin glasses with the one proposed in this paper is needed in order to explain the apparent difference in their scaling of time to solution exhibited by our numerical results.

We target the implementation of a low-power system with maximum power supply of 5W using a XCKU040 Kintex Ultrascale Xilinx FPGA integrated on an Avnet board. The implemented circuit can process Ising problems of up to at least 1100 spins fully connected and of more than 2000 spins sparsely connected within the 5W power supply limit. Data are encoded into 18 bits fixed point vectors with 1 sign, 5 integer and 12 decimal bits to optimize computation time and power consumption. An important feature of our FPGA implementation of CAC is the use of several clock frequencies to concentrate the electrical power on the circuits that are the bottleneck of computation and require high speed clock. For the realization of the matrix-vector multiplication, each element of the matrix is encoded with 2 bits precision ($w_{ij}$ is $-1$, 0 or 1). An approximation based on the combination of logic equations describing the behavior of a multiplexer allows to achieve $10^{4}$ multiplications within one clock cycle. The results of these multiplications are summed using cascading DSP and CARRY8 connected in a tree structure. Using pipelining, a matrix-vector multiplication for a squared matrix of size $N$ is computed in $2+5\frac{\text{log}(N-4)}{\text{log}(5)}+(\frac{N}{u})^{2}$ clock cycles (see Supplementary Materials S2.4) at a clock frequency of 50MHz with $u=100$ which is determined by the limitation of the number of available electronic component of the XCKU040 FPGA. The block size $u$ can be made at least 3 times larger using commercially available FPGAs, which implies that the number of clock cycles needed to compute a dot product can scale almost logarithmically for problems of size $N\approx 1000$ (see Supplementary Materials S2.4 for discussions of scalability) and that the calculation time can be further significantly decreased using a higher-end FPGA. The calculation of the nonlinearity $f_{i}$ and error terms is achieved at higher frequency (300MHz and 100Mhz) using DSP in $8+(N/u)$ and $9+(N/u)$ clock cycles, respectively. In order to minimize energy resources and maximize speed, the nonlinear and error terms are calculated multiple times during the calculation of a single matrix-vector multiplication (see Supplementary Materials S2).

Prior to computing the benchmark on the Sherrington-Kirkpatrick instances, the parameters of the system (see Supplementary Materials S2.8) are optimized automatically using Bayesian optimization and bandit-based methods[69]. The automatic tuning of parameters for some previously unseen instances is out of the scope of this work but can be achieved to some extent using machine learning techniques[70]. Sherrington-Kirkpatrick instances used in this paper are available upon request. The GSET instances are available at https://web.stanford.edu/ yyye/yyye/Gset/.

The authors thank Salvotare Mandra for providing results of the PT algorithm. This research is partially supported by the Brain-Morphic AI to Resolve Social Issues (BMAI) project (NEC corporation), AMED under Grant Number JP20dm0307009, UTokyo Center for Integrative Science of Human Behavior(CiSHuB), and Japan Science and Technology Agency Moonshot R&D Grant Number JPMJMS2021

Parameter values used for solving instances of the GSET are shown in Tab. 2.