Taming a nonconvex landscape with dynamical long-range order: Memcomputing Ising benchmarks

Recent work on quantum annealing has emphasized the role of collective behavior in solving optimization problems. By enabling transitions of clusters of variables, such solvers are able to navigate their state space and locate solutions more efficiently despite having only local connections between elements. However, collective behavior is not exclusive to quantum annealers, and classical solvers that display collective dynamics should also possess an advantage in navigating a nonconvex landscape. Here we give evidence that a benchmark derived from quantum annealing studies is solvable in polynomial time using digital memcomputing machines, which utilize a collection of dynamical components with memory to represent the structure of the underlying optimization problem. To illustrate the role of memory and clarify the structure of these solvers we propose a simple model of these machines that demonstrates the emergence of long-range order. This model, when applied to finding the ground state of the Ising frustrated-loop benchmarks, undergoes a transient phase of avalanches which can span the entire lattice and demonstrates a connection between long-range behavior and their probability of success. These results establish the advantages of computational approaches based on collective dynamics of continuous dynamical systems.

Figure 1

A schematic representation of instance creation. Separate frustrated loops (blue and red curves) are generated by random walking around the lattice until the walk crosses itself. Each loop has its own Hamiltonian consisting of $J_{ij}=1$ for all bonds except one with $J_{kl}=-1$ such that the ground state of the loop will have one unsatisfied bond. When the loops are combined, overlapping bonds (shown in black) have a coupling $J_{ij}$ which is the sum of the contributions from each separate loop.

Figure 2

Scaling of floating point operations necessary for different solvers to reach the ground state of the 3D frustrated-loop instances as a function of the total number of spins $N$. Solution times have been converted to estimated number of floating point operations. The sequential memcomputing solver implemented in MATLAB is dubbed Falcon and was run on an Intel Xeon 6148. Varying numbers of instances were run at each size and solver depending on required computation time (see the Supplemental Material for details [39]). Comparisons with simulated annealing (SA), parallel tempering (PT), and IBM CPlex run on an Intel Xeon E5430 are also shown. All calculations were performed on a single core. The solid lines are the best fits of the 95th quantile time to solution for all four solvers. The exponential fits have the following parameters: for IBM CPlex, $b=0.12$ and $c=0.46$, for SA, $b=0.069$ and $c=0.67$, and for PT $b=0.32$ and $c=0.46$.

Figure 3

An arbitrary constraint (C) satisfaction problem expressed as a factor graph can be translated into an electrical circuit with memory by considering the effect of each constraint on the site $i$. $v_R$ is a voltage generator, and $g_R$ is the conductance of a standard resistor. $v_M$ is a voltage generator, and $g_M$ is the conductance of a resistor with memory.

Figure 4

During the dynamics, memristive variables $x_{ij}$ act to interpolate between gradient dynamics on (a) the original interaction, when a constraint is unsatisfied $(x_{ij}\to 1)$, and (b) the “rigid” interaction when the constraint is satisfied $x_{ij}\to 0$. This removes the energy barrier between the two states that satisfy the interaction, allowing transitions of large clusters of variables. The satisfied configurations are denoted for the ferromagnetic case, $v_i=v_j$, and below each surface we give the effective energy function on which the voltages follow gradient dynamics.

Figure 5

When the model of Eqs. (11) and (12) is simulated for a two-dimensional instance ($L=15$) under a separation of timescales ($\beta =\frac{1}{400}$), the voltage trajectories (a) evolve under a series of sharply defined avalanches due to the slow motion of the memristors (b) modifying the clause weights. In panel (c) we have plotted the energy (left axis) without the influence of memristors (red, $\beta =0$) and with them (black, $\beta =\frac{1}{400}, \gamma =0.65$) showing that the motion of the memory variables allows the system to reach a far lower energy, and ultimately the ground state. The sizes of the avalanches (c, right axis) are plotted as gray bars, showing that their size grows over the course of the simulation until a large avalanche brings the system to its ground state. The avalanches are depicted in the inset in red with the rightmost inset corresponding to the largest avalanche in the run.

Figure 6

Spatial correlations, $\langle \mathrm{\Delta }\rangle \left(r\right)$, among voltages or spins calculated from the orthant dynamics in the slow driving limit of model (11) and (12) for the largest avalanches in two and three dimensions and for different lattice sizes. The correlations take a finite value all the way to the lattice edge, indicating that the largest avalanches are extensive. As the system size increases the values appear to saturate to a dimension-dependent value for this instance class.

Figure 7

Correlations and success probability (inset) as the minimum memristive value $x_{\min }$ is varied. When $x_{\min }=0$, no barrier remains between the two satisfied states of the interaction, and the avalanche configurations are most long-ranged. As $x_{\min }$ increases, both the correlations and the success probability decay rapidly such that by $x_{\min }=0.3$ no instances are solved.

Figure 8

Correlations and success probability (inset) as the magnitude of the rigidity terms are varied. In contrast to varying $x_{\min }$ (as in Fig. 7), which introduces a barrier between the two satisfied states, adjusting the magnitude of $R_{\mathrm{Lim}}$ maintains the ability of interacting variables to freely transition between their satisfied states. We observe a wide range in which $R_{\mathrm{Lim}}$ has a negligible effect on both the magnitude of avalanche correlations and the success probability.