Efficient combinatorial optimization by simulating quantum derived parallel annealing in memristor crossbar

Abstract

A quantum derived parallel annealing method is disclosed that achieves full degree of parallelism and improves understanding quality, resulting in significant speed and energy improvements when implemented in analog memristor crossbars. Tasks are experimentally addressed, including unweighted and weighted maximum cutting and traveler issues using integrated memristor chips. The claimed method effectively utilizes the natural parallelism provided by the memristor technology, simulated conductance states and many-to-many connections, and thus exhibits significant aging and energy efficiency improvements compared to previous simulated annealing and Isinin machines implemented in other technologies, such as for example, as well as for use in the field of memristor technology, as well as in the field of memristor technology, as well as in the field of memristor technology, as well as in the field of memristor technology. The method has a huge potential for solving a complex optimization problem with higher efficiency.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to U.S. Provisional Patent Application Serial No. 63 / 508,260 filed on June 14, 2023, and the disclosure of the U.S. Provisional Patent Application is incorporated by reference in its entirety. Technical Field

[0002] The present invention relates generally to the field of combinatorial optimization. More specifically, the present invention relates to providing an alternative approach to solving combinatorial optimization problems with higher efficiency by simulating quantum-derived parallel annealing in a memristor crossbar switch. Background Art

[0003] Combinatorial optimization problems (COPs) are ubiquitous in social life and industry, with various applications in fields including but not limited to computer science, engineering, chemistry, logistics, and economics.

[0004] Unfortunately, obtaining an exact solution is very challenging due to the combinatorial explosion that occurs as the problem size increases. Therefore, specialized hardware such as Ising machines are increasingly sought after to efficiently solve the COP by mapping the COP to the Ising model, which is a statistical model that describes a physical system consisting of interacting spins that tend to evolve to the lowest system energy.

[0005] Although classical computers can simulate the process, their efficiency is significantly limited due to their digital and serial nature and the characteristics of separate memory and processing units.

[0006] To overcome these limitations, researchers have explored a variety of analog technologies, including superconducting qubits, coherent light, CMOS oscillators, nano-oscillators, and memristor-based analog Ising machines. Among them, memristor-based machines offer particular promise due to their speed and energy efficiency, natural many-to-many connectivity, and compatibility with the electronic computing ecosystem.

[0007] Nevertheless, existing memristor-based demonstrations have not fully exploited the massive parallelism and analog storage / processing features of memristor crossbar arrays. Previous work has relied primarily on simulated annealing and its variants to obtain optimal solutions, and is therefore limited by the serial update nature of simulated annealing. Therefore, in these implementations, the memristor crossbar performs only one vector-vector dot product operation at a time, rather than a vector-matrix multiplication, resulting in the large amount of natural parallelism provided by the memristor crossbar array.

[0008] Adiabatic annealing, another annealing method with a solid theoretical background, has been successfully implemented in quantum systems and has recently demonstrated its effectiveness with a classical memristor-based platform, thus shedding light on possible parallel updating.

[0009] However, in this pioneering work, simulated annealing is still required in addition to adiabatic annealing during the solution process, and thus it does not escape its serial update limitation. In addition, most previous works only present binarized memristor conductances, rather than analog values, to represent arbitrary coupling strengths when representing problems, thus limiting the complexity of the problems that can be solved. Summary of the invention

[0010] In view of the existing challenges in the art, an object of the present invention is to provide a solution to solve combinatorial optimization problems with higher efficiency by simulating quantum adiabatic annealing using classical analog memristor crossbar switches.

[0011] According to a first aspect of the present invention, there is provided herein a quantum-derived parallel annealing processor for analyzing and solving combinatorial optimization problems with improved efficiency, comprising: Memristor crossbar switch array; Driver; multiplexer; Transimpedance amplifier; analog-to-digital converters; and Digital logic; The crossbar switch array may be formed in a single transistor single memristor configuration or in an alternative configuration.

[0012] According to one embodiment of the present invention, solving a combinatorial optimization problem includes: representing spin configurations as discrete values; using simulated variables to characterize intermediate spin states; and implementing a time-dependent system Hamiltonian; and updating the simulated intermediate spin variables using a through-the-wire estimator algorithm based on the gradient of the system Hamiltonian calculated from the actual spin configuration.

[0013] According to another embodiment of the present invention, the Ising coupling between each spin pair is stored as an analog conductance value in a memristor crossbar switch to achieve a single-step in-situ calculation of the gradient of the Ising Hamiltonian. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] For a better understanding of the various embodiments described herein, and to more clearly show how these various embodiments may be implemented, reference will be made, by way of example, to the accompanying drawings which illustrate at least one example embodiment and which will now be described.

[0015] Figure 1 The key characteristics of quantum-derived parallel annealing and its differences from simulated thermal annealing are demonstrated. Simulated thermal annealing (also called simulated annealing) has a serial update property and is different from the quantum-derived parallel annealing of the present invention in terms of annealing strategy, update scheme, and hardware implementation method.

[0016] FIG. 2A to FIG. 2F An experimental demonstration of QPA in a memristor crossbar array is shown. Figure 2A A schematic diagram of the experimental setup is provided. A general-purpose computer controls the integrated memristor chip that stores the Ising coupling matrix values ​​and computes the gradient in the analog domain, and implements part of the spin update in the digital domain. On the right is an optical image of the integrated memristor chip, which consists of multiple 64×64 single-transistor single-memristor (1T1M) arrays and necessary peripheral circuits including drivers, multiplexers (MUXs), transimpedance amplifiers (TIAs), and analog-to-digital converters (ADCs). Figure 2B The Ising coupling matrix of the target 64-node 50% density Max-Cut problem is shown, and Figure 2C The experimental readout of conductance values ​​after programming Ising coupling onto a memristor crossbar array is shown. “0” is programmed to target 0 μS, while “-1” is programmed to target 150 μS. Figure 2D and Figure 2E The success probability of QPA, SA and DHNN is ( Figure 2D ) and the average solution gap across different total iterations (solution gap = optimal maximum cut - obtained cut) ( Figure 2E ). One iteration represents one attempt to update the spin. The blue dashed line shows the simulation results of the QPA, which show strong agreement with the experimental data. Each data point is calculated based on 100 different trials based on random initial states. Figure 2F The relationship between the probability of success and the problem density for different techniques including QPA, SA, CIM and D-wave of the present invention is shown. The shaded area in the plot shows the interquartile range (25%-75%). For each problem density, 20 randomly generated problem examples are solved. 1000 trials are performed for each problem example to calculate the probability of success. For QPA and SA, the annealing time is 1000 iterations, which takes 1,000 μs if it is assumed that each iteration takes 1 μs. For CIM and D-wave, the annealing time is also 1,000 μs.

[0017] 3A to 3H A weighted maximum cut problem with simulated conductance values ​​for Ising coupling is demonstrated. Figure 3A The normalized Ising coupling matrix (-J / max(-J)) for a randomly generated 64-node many-to-many link maximum cut problem is shown. Each link strength is randomly selected from a 16-bit integer (0 to 65,535). Figure 3B The experimental readout of the conductance matrix after programming Ising coupling into a memristor crossbar array is demonstrated, where the maximum coupling strength is programmed to target 150 μS. Figure 3CThe distribution of programming error (difference between experimental readout conductance and target conductance) is shown, where the dashed line illustrates a Gaussian distribution fit with mean 0.29 μS and standard deviation 2.36 μS. Figure 3D The distribution of the calculated error (difference between experimental result and expected value) of the simulated memristor crossbar switch is shown, where the dashed line indicates a Gaussian distribution with mean 0.01 and standard deviation 0.57. 1000 random input vectors were generated for vector-matrix multiplication, where each entry was randomly taken from +1 or -1 to obtain the distribution. The expected value was calculated by software with full floating point precision. Figure 3E The calculated cut values ​​during the solution process using QPA, DHNN and SA schemes are shown. The goal is to obtain the maximum possible cut value. The cut values ​​shown here are normalized by dividing the obtained cut value by the optimal maximum cut. The light-colored line represents 100 different trials according to random initial spin configurations, and the dark-colored line represents the average cut value. The annealing time for all three methods is 1000 total iterations. Figure 3F The corresponding solution cut distributions of the three methods are shown. Figure 3G The effect of the change in conductance of the memristor is shown. For each data point, 20 randomly generated problem examples with connection strengths randomly taken from 16-bit integers are solved 1000 times with random initial spin configurations. The blue lines with square, circle, and pentagon markers represent the success probabilities of the optimal maximum cut, 99.9% maximum cut, and 99.5% maximum cut, respectively. Figure 3H The scaling trend of the time to solution (TTS) with the problem size is shown. For each problem size, 20 randomly generated problem examples with connection strengths randomly taken from 16-bit integers are solved 1000 times with random initial spin configurations. The shaded area in the plot shows the interquartile range, the square markers indicate the median, and the dashed line shows the The fitting of

[0018] FIG. 4A to FIG. 4E The traveling salesman problem is shown. Figure 4A The target N-city traveling salesman problem is presented, with solid circles marking the locations of Hong Kong, London, Johannesburg, Tokyo, Dubai, Sao Paulo, Sydney, Singapore, and New York. The goal is to find a route that traverses all cities once with the minimum distance. Figure 4B A representation of the visit order in terms of a matrix (N-1)×(N-1) is shown, where the rows represent cities and the columns represent the visit order. Each row and column can have only one element set to one, representing a valid visit order. City 1 (in this case Hong Kong) is always selected as the first city to be visited. The matrix representation in the plot is Figure 4A The example valid visit order shown in is Hong Kong>Singapore>Dubai>London>New York>Sao Paulo>Johannesburg>Sydney>Tokyo>Hong Kong. Figure 4CThe normalized Ising coupling matrix for the target problem (top) and the experimental readout conductance values ​​after programming the Ising coupling into the memristor crossbar array (bottom) are shown, respectively. The maximum coupling strength is programmed to target 150 μS. Figure 4D The evolution of the Ising Hamiltonian during the solution process is shown. The dark line represents the average cut value for 100 different trials (light lines) according to random initial states. 1000 total iterations are used for all three methods. Figure 4E A comparison of the quality of solutions obtained by QPA, SA, and DHNN is presented. Three performance metrics are compared: number of valid trips, minimum distance, and average distance. The number of valid trips represents the number of solutions that comply with the constraint that each city should be visited and visited once (the larger the better). The minimum distance refers to the shortest travel distance among 100 trials, and the average distance represents the average of the travel distances of the valid trips (the smaller the better).

[0019] Figure 5 Different H with different total number of iterations are shown 初始 The average success probability of the function. The curves with different colors represent the 初始 situation, situation, situation, situation, situation and For each data point, 20 60-node 50% density unweighted maximum cut problem examples are solved. For each problem example, 100 trials with random initial states are performed to calculate the success probability. The data in the plot shows the average success probability for the 20 problem examples. It can be observed from the plot that The symmetric function of performs better than the asymmetric function. i = 0, they all have ground states that are easy to find. and The successful implementation of also means that any function can serve as the initial Hamiltonian, as long as its ground state can be easily solved. According to empirical observations, among these three functions, The performance is the best, so it is selected as H 初始 Nevertheless, there may be other well-designed functions that can enhance performance, which requires further exploration.

[0020] Fig. 6A and 6BA comparison of different update techniques in the unweighted max cut problem is shown. a Average success probability versus total number of iterations. b Average solution gap (defined as the difference between the solved cut value and the maximum cut value) versus total number of iterations. Five different update techniques are compared: the original update algorithm used in the paper (denoted as the baseline in the figure), an update technique that does not discretize the proxy spin values ​​when computing the gradient (no STE), an update technique that does not add a momentum term, an update technique without momentum and without STE, and an update technique that does not clip the proxy spin values ​​to between -1 and 1. For each data point in the plot, 20 60-node 50% density unweighted max cut problem examples are solved. For each problem example, 100 trials are performed based on random initial states. From the plots, it can be seen that all three techniques STE, momentum and clipping are crucial to the performance.

[0021] When the total number of iterations is small, the performance of the case without STE is as good as that of STE. However, in the case of larger total iterations, the performance decreases. This is because in the case without STE, during the entire solution process, no real binary spins are involved, and the simulated spin agent cannot accurately reflect the results of the binary spin, which causes a mismatch [1]. Further updates to the simulated spin agent lack the ability to reduce the energy calculated by the binary spin. For the case without momentum, the system behaves similar to the greedy gradient descent algorithm and annealing fails. For this case, local minimum solutions can be easily obtained, however, it is almost impossible to find the global minimum. The role of momentum in quantum-derived annealing is unclear, but it may be related to the characteristics of STE. Because the performance of ordinary stochastic gradient descent is much worse than those with momentum in the field of binary neural network (BiNN) training. For the case without momentum and without STE, annealing returns to normal. The convergence rate is relatively slower than that of the case without momentum. For the case without momentum and without STE, annealing works normally, but the convergence rate is relatively slower than that of the case without momentum. For the case of infinite amplitude, the simulated spin agent has no restrictions and can grow infinitely in one direction during the solution process. Then, if it needs to go in another direction later in the solution process, it takes much longer to accumulate the gradient until it can change its sign, and thus greatly reduces the solution speed. For STE, momentum and clipping, they are all common techniques used in BiNN training. Although there are more advanced optimizers such as Adam and other precision training methods that are widely used in the field of neural networks, they have the potential to be applied to combinatorial optimization, but due to the different nature of the target problem, their applicability needs to be further studied.

[0022] Fig. 7A and Figure 7B The overall experimental setup is shown. Fig. 7AAs shown in Figure 1, a PC is used to control the integrated memristor chip via a microcontroller (left), and a printed circuit board (PCB) is developed to measure the chip (middle). On the right, an optical image of the integrated memristor chip is shown, which consists of three 64×64 1T1M memristor arrays and peripheral circuits, including a transimpedance amplifier (TIA), a multiplexer (MUX), an analog-to-digital converter (ADC), and a driver. The CMOS peripheral circuits are taped out with TSMC's 180nm technology, while the memristor device is manufactured in the back-end of the upper production line with an internal clean room. Figure 7B An optical image of a memristor crossbar array is shown (left); a SEM image of a crosspoint memristor device (center); and a lateral TEM image of a 50nm×50nm Pt / Ta / TaOx / Pt memristor device (right).

[0023] Fig. 8A and Figure 8B The performance of a memristor device in integrated hardware is demonstrated. Fig. 8A A DC IV measurement of 100 set and reset cycles is shown, showing small cycle-to-cycle variation. The average is shown by the black line. Figure 8B Linearity testing of different conductance states is demonstrated, showing excellent IV linearity. Figure 8C Shown is a hold test after partitioning the entire 64×64 array into 16 subarrays, each programmed to target a conductance of 0 μS to 150 μS with a step size of 10 μS. The mean is shown by the black line, while the shaded area indicates the interquartile range.

[0024] Fig.9A and Fig. 9B The influence of the initial λ value and its relationship to the maximum eigenvalue of the Ising coupling matrix are shown. Fig.9A The Ising coupling matrix is ​​shown. The maximum eigenvalue is 7.76. Fig. 9B The success probability vs. initial λ for the problem example in (a) is shown. The black dashed line shows the location of the initial λ that is equal to the maximum eigenvalue of the Ising coupling matrix.

[0025] In the QPA of the present invention, the coefficients of the initial Hamiltonian should start with sufficiently large coefficients so that the system can start with a function (i.e., usually a convex function) that can easily find the ground state. i The calculated gradient is used to update the simulated agent spin x i By replacing the binary spin configuration σ between the Ising Hamiltonians with the simulated spin configuration x, we obtain the following:

[0026] In this case, if λI-J ≥ 0, then H系统 becomes convex and when all x i = 0. Therefore, to ensure a convex function, λ should start with a value greater than the maximum eigenvalue of the Ising coupling matrix J. However, it cannot be too large, because for a larger initial λ, more iterations will be wasted and thus performance will be affected. From the above plot, it can be seen that the best performance point is at the position where the initial λ value is greater than the maximum eigenvalue of the Ising coupling matrix.

[0027] Fig. 10A and Fig. 10B Demonstrate the effect of step length. Fig. 10A shows a plot of success probability versus step length; and Fig. 10B shows a plot of the average solution gap versus step size. For this plot, 100 trials were run for Figure 2C The dotted line shows the simulation results using the experimentally verified memristor model, which fits the experimental data well. As can be seen from the plot, when the step size is relatively small, increasing the step size can lead to improved system performance because each iteration advances further in the gradient direction. However, as the step size becomes large enough, the system will diverge. This is consistent with the experience in neural network training. In fact, more advanced adaptive changes in step size are usually implemented in neural network training. The selection of optimal hyperparameters including step size and momentum term coefficient is a huge and complex problem that requires further study in the future.

[0028] Fig.11A and Fig. 11B The evolution of intermediate values ​​during the experimental solution of the weighted maximum cut problem is shown. Fig.11A The evolution of a classical "superposition" (simulating a spin agent) is shown. Fig. 11B The evolution of momentum is shown. During the early stages of the solution process, the classical "superposition" oscillates rapidly between -1 and 1 due to the high influence of the initial Hamiltonian, and eventually converges to -1 or 1 to match the binary spin configuration.

[0029] Fig. 12A and Fig. 12B Shown is a plot of the success probability of computing TTS for the weighted maximum cut problem. Fig. 12A Shown is a plot of the probability of success versus problem size for simulated annealing. Fig. 12B Shown is a plot of the probability of success versus problem size for quantum-derived annealing.

[0030] Fig.13A and Fig. 13B Scaling trends are shown for the unweighted maximum cut problem. Fig.13A Shown is a plot of the probability of success versus problem size for simulated annealing. Fig. 13BShown is a plot of the probability of success versus problem size for quantum-derived annealing. Fig. 13C A plot of solution time versus problem size is shown. For each problem size, 20 randomly generated problem examples were solved 1000 times with random initial spin configurations. The shaded area in the plot shows the interquartile range. The square marker indicates the median.

[0031] FIG. 14A to FIG. 14H A TSP example and travel path comparison are shown. Fig.14A An example of the target problem solved in this paper is shown. Fig. 14B The corresponding distance matrix between each city in km is shown. The distance between each two cities is the shortest distance along the surface of the earth and calculated by the coordinates of the city centers. Fig. 14C shows a spin-value matrix representation of the minimum distance travel path found by the QPA, and Fig.14D Its corresponding travel path in the map is shown. Fig.14E shows a spin-value matrix representation of the minimum distance travel path found by SA, and Fig.14F Its corresponding travel path in the map is shown. Figure 14G shows a spin-value matrix representation of the minimum distance travel path found by the DHNN, and Fig.14H Its corresponding travel path in the map is shown.

[0032] FIG. 15A to FIG. 15C Simulation results for a larger problem solved in TSPLIB are shown. FIG. 15A to FIG. 15C Results for solving the burma14 problem, the ulysses16 problem, and the gr17 problem are shown, respectively. The best in the plot represents the best known solution provided by the library. For each problem, 5000 iterations were performed for both QPA and SA by 100 trials based on random initial states. DETAILED DESCRIPTION

[0033] In the following description, the method and modification procedure of the aptamer as a target detector are described as a preferred example. It should be understood by those skilled in the art that modifications, including additions and / or substitutions, can be made without departing from the scope and spirit of the invention. Certain details may be omitted so as not to obscure the invention; however, this disclosure is written to enable those skilled in the art to practice the teachings herein without undue experimentation.

[0034] Various processes will be described below to provide examples of at least one embodiment of the claimed subject matter. The embodiments described below do not limit any claimed subject matter, and any claimed subject matter may cover processes or systems different from those described below. The claimed subject matter is not limited to processes or systems with all the features of any process or system described below, or is not limited to the features common to multiple or all processes or systems described below. The process or system described below may not be an embodiment of any claimed subject matter. Any subject matter disclosed in the process or system described below that is not claimed in this document may be the subject matter of another protective instrument, for example, a continuing patent application, and the applicant, inventor, or owner does not intend to abandon, deny, or dedicate any such subject matter to the public through the disclosure in this document.

[0035] In addition, it should be understood that for simplicity and clarity of illustration, reference numerals may be repeated between the figures to indicate corresponding or similar elements where deemed appropriate. In addition, numerous specific details are set forth to provide a thorough understanding of the embodiments described herein. However, it should be understood by those of ordinary skill in the art that the embodiments described herein may be practiced without these specific details. In other examples, well-known methods, procedures, and components are not described in detail to avoid confusing the embodiments described herein. In addition, the description is not to be construed as limiting the scope of the embodiments described herein.

[0036] Background on quantum annealing in quantum-derived processors

[0037] Quantum annealing is a simulation computing technique for solving optimization problems by exploiting quantum phenomena. Quantum annealing uses the principles of quantum mechanics to explore the solution space of a given optimization problem more efficiently than classical methods.

[0038] In the context of quantum annealing, spin refers to a quantum mechanical entity that can represent a binary variable. Each spin can be in one of two states, which are usually labeled "speeded up" or "slowed down," which can be analogous to the values ​​0 and 1 in classical computing.

[0039] The Ising model is a mathematical model for describing the behavior of spins interacting with each other in a system. The Ising model is often used in statistical mechanisms and condensed matter physics to study phenomena involving phase transitions. In the context of quantum annealing, the Ising model is suitable for representing optimization problems; and the interactions between spins in the Ising model are called Ising couplings. These interactions are represented by coupling coefficients, which determine the strength with which spins affect each other. In optimization problems, these coupling coefficients encode the relationship between the variables that need to be optimized.

[0040] In contrast to the description of the total energy of a quantum system, in quantum annealing, the Hamiltonian is formulated to encode the optimization problem as an Ising model with certain constraints. It consists of two parts: a "cost function" term, which represents the objective function of the optimization problem; and a "tuning" term, which controls the annealing schedule, which refers to the evolution of the system over time. Initially, the system is prepared in a state where the tuning term dominates, allowing the system to explore the solution space extensively. Over time, the tuning term gradually decreases, while the cost function term becomes more dominant, guiding the system towards the optimal solution.

[0041] Throughout the annealing process, measurements are performed on the system to extract information about its state. These measurements provide probabilistic information about the solution to the optimization problem encoded in the system's quantum state.

[0042] At the end of the annealing process, the system ideally enters a state that represents the solution to the optimization problem. Measurements can be performed on this final state to extract the solution with a certain probability.

[0043] Memristor crossbars and arrays can be used as hardware implementations for quantum annealing systems. While they are not inherently quantum devices themselves, they can be used as part of classical control and readout systems in quantum annealing settings.

[0044] A memristor is a two-terminal electronic device whose resistance depends on the history of charge that has flowed through it. This unique property allows memristors to store and process information in a non-volatile manner, making memristors a promising candidate for future memory and computing technologies.

[0045] Furthermore, in the context of memristors, a crossbar architecture refers to a grid-like arrangement in which rows and columns of memristors intersect. Each intersection represents a memory cell or a computational cell. By applying appropriate voltages to the rows and columns, it is possible to control the state of individual memristors and perform operations such as matrix-vector multiplication, which is fundamental to many computational tasks including optimization problems.

[0046] Memristor arrays are large collections of memristor crossbars arranged in a grid-like fashion. These arrays can be scaled up to thousands or even millions of memristors, allowing for massively parallel computation and storage. In the context of quantum annealing, memristor arrays can store and manipulate the coefficients of an Ising model representing the optimization problem being solved.

[0047] While the Ising model coefficients encode the relationship between the variables in the optimization problem, memristor arrays can be used to efficiently store these coefficients. The Ising model can be mapped onto a physical memristor array by setting the resistance of each memristor in the array to a value corresponding to the coefficient it represents.

[0048] The memristor crossbar can perform the computational operations required for quantum annealing, such as matrix-vector multiplication. By applying appropriate voltages to the rows and columns of the crossbar, it is possible to perform a dot product between the Ising model coefficients and the current state of the quantum annealing system. This allows efficient calculation of the energy landscape of the optimization problem, which is crucial for guiding the annealing process.

[0049] Memristors can also be used to provide feedback and control signals to quantum annealing systems. For example, they can be used to read out the current state of a quantum system or adjust an annealing schedule based on measured performance.

[0050] According to a first aspect of the present invention, there is provided a quantum-derived parallel annealing processor for analyzing and solving combinatorial optimization problems with improved efficiency, comprising: an integrated memristor chip including a plurality of 64×64 single-transistor single-memristor crossbar switch arrays; a driver; a multiplexer; a transimpedance amplifier; and an analog-to-digital converter; and a general-purpose computer as a control device for the integrated memristor chip.

[0051] In this configuration, when a memory operation begins, the driver may apply appropriate voltages or currents to the rows and columns of the 64×64 single-transistor single-memristor crossbar switch array for reading and writing data. Correspondingly, the multiplexer is configured to select and connect specific rows or columns from multiple rows and columns to read or write data. Subsequently, the transimpedance amplifier converts the current signal from the 64×64 single-transistor single-memristor crossbar switch array into a voltage signal for subsequent signal processing. In addition, the analog-to-digital converter is configured to convert the amplified analog voltage signal into a digital signal for digital processing and storage. In one embodiment, the quantum-derived parallel annealing processor further includes a digital logic component configured to perform digital signal processing, control the operation of the entire system, and handle data reading, writing, and processing.

[0052] The above example demonstrates how a 64×64 single-transistor single-memristor crossbar array can be combined with other components to implement a complete memory system. This system can be used for memory storage and other applications that require efficient data storage and processing based on the principles of the Ising Hamiltonian and spin states, and in the example of the present invention, for solving and analyzing combinatorial optimization problems.

[0053] According to one embodiment of the present invention, solving a combinatorial optimization problem includes: representing a spin configuration as a discrete value of 1 or -1; using simulated variables to depict intermediate spin states equivalent to a classical "superposition" of spins; and implementing a time-dependent system Hamiltonian that starts with an initial Hamiltonian dominating and gradually transitions toward an Ising Hamiltonian; and using a through-the-wire estimator algorithm to update the simulated intermediate spin variables based on the gradient of the system Hamiltonian calculated based on the actual spin configuration.

[0054] According to another embodiment of the present invention, the Ising coupling between each spin pair is stored as an analog conductance value in a memristor crossbar switch to achieve a single-step in-situ calculation of the gradient of the Ising Hamiltonian.

[0055] Quantum derived parallel annealing (QPA)

[0056] The combinatorial optimization problem (COP) can be reformulated as the Ising model, where the spin configuration evolves towards the ground state of a certain energy function called the Ising Hamiltonian, expressed as where σ={σ1,σ2,…,σ i ,…,σ N} is the spin configuration vector encoding the solution to the problem, where each component σ i is the binary value {-1,1} representing the two states of spin-up and spin-down. J is a square and symmetric coupling matrix with size N×N, where each element represents a ferromagnetic or antiferromagnetic connection between two spins. h is a local field term introduced for generality.

[0057] In physical quantum systems, the annealing process for solving such models is achieved by adiabatic evolution. It is based on a simple Hamiltonian H whose ground state can be easily found. 初始 (For example, the transverse field Hamiltonian ) and gradually shifts to the Ising Hamiltonian Expressed as H 系统 (t) = A(t)H 初始 +B(t)H 伊辛 (2)

[0058] Where A(t) gradually changes from 1 to 0, and B(t) gradually changes from 0 to 1. Ideally, the physical system can always maintain the minimum energy state, and thus the system can eventually evolve to the ground state of the Ising Hamiltonian and solve the corresponding COP.

[0059] However, in this quantum version of adiabatic annealing, the spins are represented by matrices rather than discrete values, which are not suitable for simulation in this current classical analog memristor crossbar. Therefore, it is adapted for easier implementation in this analog memristor crossbar, and a new classical version of adiabatic annealing is proposed, in which the spins σ i Represented by discrete values ​​1 or -1. In order to perform annealing, a simulation variable x is introduced i to represent the intermediate spin states, and the simulated variables can be considered as a classical "superposition" of spins. The real spin configuration σ i Treated as symbol x i .

[0060] Similar to the quantum version, the time-dependent system Hamiltonian is implemented as follows: H 系统 (t) = λ(t)H 初始 +H 伊辛 (3) where λ(t) is the time-dependent coefficient, which starts with a sufficiently large value and gradually decreases to 0 during the solution process. 初始 is the initial Hamiltonian and can be any function whose ground state (global minimum) can be easily found. In this work, we choose It is one of the simplest convex functions (H 初始 The different options and their comparison can be found in Figure 5 In the solution process, the system is firstly 初始 dominant, and will gradually shift to H 伊辛 , thus generating H 伊辛 However, unlike real physical systems that naturally tend to evolve towards the lowest energy state during adiabatic evolution, adiabatic Hamiltonian shifts require the development of algorithms to achieve manual updating of the spins and thus lowering of the energy.

[0061] To address this challenge, gradient descent, a popular optimization technique used in artificial neural network training, is applied to help the system dynamically evolve to a lower system Hamiltonian during adiabatic shift. However, in this case, ordinary gradient descent cannot be applied due to the discrete nature of the spin configuration. Each update can only flip the sign of the spin configuration, resulting in divergence. To overcome this limitation and derived from the training of binary neural networks (BiNNs), the straight-through estimator (STE) algorithm is adopted, which is a modified version of gradient descent that has shown great success in the field of neural networks. The key idea behind STE is to introduce full-precision "latent" weights as a proxy for binary weights, which are binarized in the forward and backward paths to calculate the gradient value. This gradient value is then directly used as an estimator to update the full-precision "latent" weight. Interestingly, the analog intermediate spin variable x introduced earlier to represent the classical "superposition" of spin states i Shares similar properties as the "latent" weight, and can therefore serve as a proxy for spin.

[0062] To implement the STE algorithm in an Ising machine, the simulated spin variable x is first projected onto the binary domain using a symbolic function to obtain the real spin configuration σ. The real spin configuration is then used to compute the gradient of the system Hamiltonian. Specifically, the gradient is given by:

[0063] This gradient is then used to update the simulated spin agent x as x(t+1)=x(t)-η*gradient (5), Where η is the step size. To improve the performance of the algorithm, two common techniques from stochastic gradient descent in neural network training have been implemented: clipping and momentum, which are widely used in stochastic gradient descent in neural network training. More techniques from modern neural networks can potentially be applied in the future to further improve performance.

[0064] Notably, the gradient update based on is applied simultaneously with the adiabatic Hamiltonian shift to dynamically find the ground state of the current system Hamiltonian and eventually evolve to the ground state of the Ising Hamiltonian. During the solution process, the Ising coupling between each pair of spins is stored in a many-to-many manner in the memristor crossbar as analog conductance values, which can be used in situ to compute the gradient of the Ising Hamiltonian in a single step. Therefore, all outputs from the memristor crossbar can be used to update the spins, thus preserving its parallel and analog nature. Figure 1 The key characteristics of this QPA implemented in an analog memristor crossbar and its differences from previous simulated thermal annealing are outlined in.

[0065] Demonstration of QPA for the maximum cut problem in a memristor crossbar switch

[0066] In order to demonstrate the solving capability of the currently presented QPA in a memristor crossbar, maximum cut was chosen as the benchmark problem. Maximum cut is a classic and widely studied NP-hard combinatorial optimization problem that is often used to benchmark Ising machines. The choice is due to its direct mapping capability to the Ising model and important practical applications including integrated circuit routing, computer vision, and data clustering. The maximum cut problem involves partitioning all vertices V of a given graph G(V,E) into two sets in a manner that cuts the maximum number of edges connecting nodes in different sets. A detailed explanation of the mapping of this problem to the Ising model is found in the methods section. In an embodiment, the Ising coupling is stored in the memristor crossbar as an analog conductance value that is programmed using an iterative write and verification method. The memristor crossbar is then used to calculate the gradient by performing matrix multiplication in the analog domain, while a portion of the spin update is performed in the digital domain ( Figure 2A All experiments were conducted on the integrated memristor chip of the present invention, which consists of a plurality of 64×64 single-transistor single-memristor (1T1M) arrays and necessary peripheral circuits, including a driver, a multiplexer (MUX), a transimpedance amplifier (TIA), and an analog-to-digital converter (ADC) ( Figure 2A In this demonstration, a personal computer (PC) is used to control the chip and update the spin ( Figure 2A ). More details on the experimental setup and electrical characterization of the memristor device can be found in Figures 7A to 7B and Figures 8A to 8C .

[0067] First, we use a 50% density (connecting 50% of the edges) unweighted 64-node max-cut problem to ensure a fair comparison with recently reported results. Figure 2B The Ising coupling matrix mapped from the problem is shown in . Since only negative Ising coupling is required in the maximum cut problem and the conductance values ​​can only be positive, the negative Ising coupling (-J) is converted into conductance values ​​in the range of 0μS to 150μS. After configuring the conductance values ​​into the hardware, the resulting conductance matrix is ​​experimentally read out as Figure 2C The consistency between the mapped matrix and the experimentally read matrix demonstrates the good performance of this device and hardware platform. Once programmed, the problem is solved using the method described previously without changing the memristor conductance value.

[0068] Then the problem is solved by QPA. Figure 2D and Figure 2EThe performance of the QPA is compared with a naive discrete-time Hopfield neural network (DHNN) and simulated annealing (SA for short) over different annealing times (total iterations) in . The DHNN is a local search algorithm that serially updates spins towards lower energy states. The SA implementation here involves adding noise with a decreasing amplitude to the DHNN, which is equivalent to the decreasing "temperature" in the classical SA implementation, and the SA implementation is commonly used in memristor-based implementations. Throughout the iterations, the noise level (standard deviation of the Gaussian distributed noise) in the SA changes linearly from 2 to 0. Similarly, for the presented QPA, λ changes linearly from 10 to 0, with a fixed step size η=0.01 used for the gradient descent portion (further discussion on the choice of λ and step size η can be found in , respectively). Figures 9A to 9B and Figures 10A to 10B ). The results show that the DHNN without annealing has a negligible chance of success due to the local search properties, while the method of the present invention clearly gives better solutions than SA in terms of both the probability of success in reaching the maximum cut value (ground state) and the average cut value with the same total iterations. This is due to the parallelism of the proposed method, which allows the system to process and protect more information in a single iteration. In addition, the annealing speed is affected by the change in the total number of iterations in this experiment, because the noise level in SA and the coefficient λ of the initial Hamiltonian in this QPA change linearly throughout the iterations. The success probability of SA will eventually reach 100% with a sufficient number of iterations, because it has been established that if the "temperature" changes slowly enough, the convergence of SA is guaranteed. The same is true for adiabatic annealing in quantum systems, where the transfer from the initial Hamiltonian to the Ising Hamiltonian should be slow enough (i.e., follow the adiabatic evolution rule) to ensure that the system remains in the ground state. Based on the experimental results that the probability of success continues to increase with the increase in the total number of iterations, this may also be the case for this classical version of adiabatic annealing.

[0069] The performance changes with different question densities are further investigated and the results are compared with those obtained using other techniques ( Figure 2F ). Due to the lack of guaranteed convergence in the absence of annealing, DHNNs are not compared in this analysis. Due to experimental hardware resource limitations, and considering the non-idealities in the experiments, a memristor crossbar switch model was experimentally validated, which has shown good agreement with experimental data from previously described experiments ( Figure 2D and Figure 2E ). The results show that QPA can achieve a higher success probability than all competing methods including SA, coherent Ising machine (CIM) and D-wave quantum annealer. In addition, since the memristor crossbar array is dense in nature and thus has many-to-many connectivity, the proposed system can easily handle problems with different densities.

[0070] Weighted Maximum Cut Solution for Many-to-Many Connections

[0071] One major difference between this approach and other techniques used in Ising machines is the use of a memristor crossbar array to store coupling strengths, which provides a unique feature: each crosspoint can be programmed to an arbitrary conductance state, thereby enabling the representation of arbitrary coupling strengths between any spins with many-to-many connectivity. Since practical problems often require more levels of coupling strength, this feature allows a system based on the method of the present invention to solve such problems without any additional hardware cost, so as to significantly improve both speed and energy efficiency.

[0072] Previous results have demonstrated that this system can be used to solve highly dense max-cut problems. In addition to the unweighted problem, the many-to-many connected weighted max-cut problem has been experimentally solved to fully exploit the analog storage and processing capabilities of the memristor crossbar array. First, a random weighted max-cut problem is generated, where the weight of each edge is randomly assigned a random 16-bit integer (from 0 to 65,535). The problem is then mapped to the appropriate conductance range (0 μS to 150 μS) and programmed onto the memristor crossbar array ( Figure 3B ). Figure 3C and Figure 3D The programming and computational accuracy are shown in . The results are compared with SA and DHNN. Figure 3E shows the evolution of the cut value during the solution process, while Figure 3F The final solution distribution is plotted in . Fig.11A and 11B The classical “superposition” evolution process is observed in . Similar to the unweighted maximum cut problem, this method obtains significantly better solutions within a certain number of iterations.

[0073] In the case of finite annealing time, as the weights become arbitrary and more energy states are introduced, the problem becomes difficult to solve and causes the performance of SA and DHNN to be worse than when they are used to solve the unweighted problem ( Figure 3E and Figure 3F ). However, the QPA of the present invention shows very good performance: 48 out of 100 trials with random initial states eventually converge to the true ground state, while SA and DHNN did not find the optimal solution even once. Although the variation of the memristor prevents it from being fully programmed to a given conductance state, the experimental results confirm that the current programming accuracy is sufficient to solve the COP with a numerically equivalent success probability. Further simulations with different conductance variations were performed to examine the impact of programming errors of the memristor crossbar array ( Figure 3G). In case of large conductance changes, the success probability of reaching the optimal solution decreases. However, in case of relatively small conductance changes (0 μS to 5 μS), the average solution quality shows negligible degradation in both success probability (blue line) and cut value (red line), and the current experimental conductance change (about 2.36 μS, as shown in Figure 2) is Figure 3C In addition, the conductance variation of memristor devices can potentially be further reduced through new material stacks and de-noising techniques to enhance the performance of future systems. For most practical application scenarios where a good enough suboptimal solution (e.g., 99.5% of the maximum cut value) is acceptable, analog computing systems exhibit better tolerance to conductance variations, as shown in Figure 3G The nearly undegraded accuracy (blue line with pentagonal markers) is shown for conductances up to 10 μS (more than four times the current experimental value). This demonstrates that the heuristics of the Ising method are particularly well suited for analog computing hardware.

[0074] In addition, the scaling of the time to solution (TTS) of different methods with the problem size is studied ( Figure 3H ). TTS is defined as the time required to guarantee a 99% probability of success in reaching the global optimal solution, which is usually achieved by performing multiple annealing runs and selecting the best result. Therefore, TTS can be quantitatively calculated as Where T ann is the annealing time required for a single run, and P is the success probability of a single run. The success probability used to calculate TTS can be found in Fig. 12A and 12B , while the scaling trend for the unweighted maximum cut problem is shown in Fig.13A and 13B It is worth noting that for QPA the results shown here are from simulations taking into account experimentally available calculation errors, whereas for SA the results are obtained with defect-free simulations, so the actual improvement for SA is larger than the one reported here. Figure 3G The results in show that, consistent with previous studies, the time complexity of both methods scales exponentially with the square root of the problem size N, i.e., where a and b are constants. Although this QPA exhibits a similar scaling trend as SA, the scaling factor is much smaller, indicating a more significant scaling advantage. For example, in the case of a problem size of N = 120, the TTS of the proposed method is 46 times lower than that of SA. Moreover, by simply extrapolating the curves, the expected speedup increases with the problem size.

[0075] Traveling Salesman Problem

[0076] In addition to its ability to solve the maximum cut problem, due to the versatility of the Ising model, the memristor-based system of the present invention also has the potential to be used to solve other types of combinatorial optimization problems (COPs). To demonstrate this, the Traveling Salesman Problem (TSP) was chosen as another classic COP benchmark task. TSP involves finding the shortest path that visits each city once and returns to the starting city, and has various applications in scheduling and routing problems. To map TSP to the Ising model, the (N-1) 2 spins, where N represents the number of cities, and the spins are arranged as a (N-1)×(N-1) matrix, where rows represent cities and columns represent the order of visits ( Figure 4A-4B ). Each row and column should have exactly one spin in the spin-up state to satisfy the constraint that all cities must be visited once and only once. For the spin variable b i Using the binary bit formula (0 or 1), where where σ i is the spin value, which is -1 or 1, indicating the spin-slowed state and the spin-fast state, respectively. This formula is used to conveniently express the energy, and the coupling matrix J and the bias h are modified accordingly. The detailed process of mapping TSP to the Ising model is presented in the method section.

[0077] Figure 4C Shown in the Figure 4A The ideal Ising coupling strength and experimental readout conductance matrix after the target problem in is mapped into conductance and configured into memristor hardware (more details about the target problem can be found in Figures 14A to 14H ). Similar to the maximum cut problem, in Figure 4D The evolution of the Ising Hamiltonian using QPA is compared with DHNN and SA. Figure 4E The final solutions obtained by different solvers after 1,000 iterations are compared in . QPA achieves a much higher solution quality with the same total number of iterations, generating more efficient runs and finding runs with smaller minimum and average distances. This demonstrates the effectiveness of this approach in solving the TSP. Further scaled simulation results on real problems solved in TSPLIB can be found in Figures 15A to 15C In addition, the benefit of the analog characteristics of the memristor device of the present invention enables a significant reduction in the number of device units required to solve the TSP. For example, to solve the 10-city TSP, the array size can be reduced from 526×200 to 81×81G.

[0078] It is worth noting that mapping the method from TSP to the Ising model may face scalability issues because the number of spins required increases with the square of the number of cities to be visited. Through advanced clustering techniques, the problem can be alleviated by breaking down the large problem into several levels of smaller problems. In this case, the speed of solution and the quality of the solution of the small problems can be crucial to the overall large problem. Therefore, the analog nature of the device combined with the parallel quantum-derived annealing of the memristor-based system of the present invention to achieve better solution quality is very suitable for such techniques. In addition, the mixed signal processing of the current embodiment becomes attractive because it is more compatible with the higher-level processing that is inevitable for mapping and clustering.

[0079] Table 1 below compares the key characteristics and performance metrics of various Ising machines for solving the 100-node dense maximum cut problem. Key characteristics include representation methods for spins and couplings, connectivity and accuracy of couplings, update and annealing mechanisms, providing a basic understanding of each technology. Performance metrics include annealing time, solution time, power dissipation, and energy efficiency. Seven different technologies are compared, including the present invention's memristor-based QPA, memristor-based Hopfield network (mem-HNN), continuous-time dynamic system based on phase-transition nano-oscillator (PTNO), Ising system based on CMOS ring oscillator (ROSC), analog bifurcation machine (SBM) running discrete analog bifurcation (dSB) on field programmable gate array to represent the current state-of-the-art digital solver, coherent Ising machine (CIM), D-wave 2000Q quantum annealer. Details on the estimated breakdown of the QPA are shown in Table 2 below. Table 1: Benchmark comparison of QPA implemented in a memristor crossbar with other state-of-the-art Ising machines for a 100-node 50% density unweighted max-cut problem Table 2: Energy and area estimation breakdown Modules <![CDATA[Area / μm 2 > Delay / ns Energy / pJ Array 137.21 <0.3 1.54 Drive 275.56 0.1 0.34 S&H 3.16 - 0.02 MUX 101.80 0.008 0.20 ADC 14250 - 2.50 sum 14767.73 - 4.60

[0080] To ensure a fair comparison with mem-HNN, the estimation is based on the 16nm technology node. For the analog-to-digital converter (ADC), an 8-bit state-of-the-art successive approximation register (SAR) ADC design was used for the estimation. For the QPA implementation, the solution time was benchmarked using the same operating frequency as the mem-HNN at 500Mhz (2ns per iteration). The power dissipation of this system was estimated by dividing the energy consumption per iteration by the time required per iteration (4.60pJ / 2ns) and the result was 2.3mW.

[0081] The performance of solving the 100-node 50% density unweighted maximum cut problem was compared between the inventive QPA implemented on a memristor crossbar array and six other Ising machines based on different technologies. For the inventive memristor-based QPA, the annealing time was set to 200 iterations (400ns), which yielded a 16% success probability. 27 runs were required to reach a 99% success probability, resulting in a solution time of 10.8μs. The solution energy was counted by multiplying the energy required for a single iteration (4.60pJ) by the total iterations required to reach a 99% success probability (5400), resulting in an energy consumption of 24.8nJ. For the memristor-based Hopfield network (mem-HNN), a hybrid update scheme was utilized, which updated 10 nodes at a single iteration. In order to achieve the optimal solution time, the annealing time was chosen to be 50 update cycles, where one update cycle is equivalent to N / 10 iterations (N is the problem size), resulting in an annealing time of 1μs. In addition, this hybrid update method breaks the guaranteed convergence condition of simulated annealing and may lead to divergence, especially for denser problems with large connectivity between nodes. More research is needed on the general applicability of this hybrid method. For systems based on phase transition nanooscillators (PTNOs), the fully analog characteristics achieve low energy consumption. However, the experimental demonstration is limited to a problem size of N=8, which raises doubts about the effectiveness of the system on larger problems, especially considering the larger parasitic resistance and capacitance for larger problems. For systems based on ring oscillators (ROSCs), although the annealing time is ultrafast compared to other techniques, the lack of its annealing procedure results in a poor success probability and therefore a relatively long solution time. Since the success probability of ROSC is not reported, an upper bound of 1% is used for estimation based on previous experience with local search algorithms. This results in 459 runs required to reach a 99% success probability and a corresponding solution time of 23μs, without considering the overhead introduced by comparing such a number of solutions. In addition, since the ROSC system only supports sparse connections, the embedding technique requires more spins to solve dense problems, which may further degrade performance. For the simulated bifurcated machine (SBM), the reported solution time is used to solve the many-to-many connection problem with both positive and negative coupling. However, since the difficulty of the problem has not changed fundamentally, it is used here for estimated comparison purposes. In addition, the energy consumption of the field programmable gate array (FPGA) is estimated based on its power dissipation, which is considered to be 200W in this estimate. For the coherent Ising machine (CIM), since the energy consumption is not reported, the estimated lower bound is estimated by considering only the energy consumption of the FPGA part (assuming the power dissipation is 200W) without considering the energy consumption of the optical part. For the D-wave 2000Q quantum system, the energy consumption is limited by cryogenic cooling, and the speed of solving dense problems is limited by spare connections.

[0082] Therefore, as clearly shown in Tables 1 and 2, by implementing the present QPA in a memristor crossbar array, the computational speed for operation is increased as the computer power consumption and dissipation during operation is reduced.

[0083] The 100-node unweighted max-cut problem was chosen as a benchmark for comparison purposes because it is commonly used in other reports. The QPA implemented on this memristor-based system achieved a solution time of 10.8 μs, 2.3 times faster than the previous state-of-the-art solver, and achieved a 4.10×10 7 The energy efficiency of solutions / second / Watt is 3.1 times greater than the previous state-of-the-art solver. This advantage is mainly due to the novel quantum-derived annealing scheme, which further exploits the parallel, many-to-many connectivity and analog nature of the memristor crossbar array. It is important to note that the 100-node maximum cut is not a limitation of the memristor-based system, as the current state-of-the-art memristor-based in-memory computation macro has 1024×512 devices in a single bank. In the case of larger arrays, the advantages brought by synchronous updates can also be amplified due to the exploitation of higher parallelism.

[0084] Mixed-signal memristor-based methods (including the memristor-based QPA and mem-HNN of the present invention) store and compute coupling terms in the analog domain while implementing spin updates in the digital domain. In contrast, CIM updates spins in the analog domain and computes coupling terms in the digital domain. Considering that the utilization of quantum mechanics demonstrated by current CIM is still unclear and can be accurately described by classical dynamics, it is believed that it may be more efficient to implement the coupling terms in the analog domain at the current stage before CIM shows quantum advantage. This is because the computational complexity of spin updates is typically O(N), where it is possible to achieve a time complexity of O(1) if custom parallel digital logic is implemented. On the other hand, the coupling computation implemented by CIM in the digital domain is a VMM operation with a computational complexity of O(N). 2 ), making it computationally significantly more intensive than spin updates.

[0085] Examples Limits and Momentum

[0086] During the solution process, after each update iteration, the simulated agent spin x is clipped between -1 and 1 using the following equation: Limit(x,-1,1)=max(-1,min(1,x)) (6)

[0087] This clipping is a common practice in Binary Neural Network (BiNN) training to prevent parameter values ​​from growing indefinitely, so that a slight change in the value will not have any effect on the result of the binarized parameter.

[0088] In order to further increase the convergence speed, momentum gradient descent is used: m(t+1)=β*m(t)-η*gradient (7) x(t+1)=x(t)+m(t+1) (8) Where β is the momentum constant, which is set to 0.99 in this paper for simplicity. The momentum m is clipped between -1 and +1 to prevent explosion. Mapping the maximum cut to the Ising model

[0089] The maximum cut problem aims to partition all nodes into two subsets. To express the problem mathematically, if the i-th node is in one subset, then use s i = 1, and if the i-th node is in another subset, use s i = -1. A is an adjacency matrix, where if there is no edge between the i-th node and the j-th node, the element is defined as A ij = 0, and if there is an edge, the element is defined as A ij =1.

[0090] The number of cuts can then be expressed as

[0091] Because i＜j A i,j is a constant defined for the problem, so maximizing the number of cuts is equivalent to minimizing -∑ i＜j (-A i,j )s i s j , and thus the problem can be mapped to the Ising Hamiltonian by setting J = -A, without the need for the local field term h. Mapping TSP to the Ising model

[0092] To map the N-city TSP problem, N spins are needed. Each spin is represented by a binary bit (0 or 1) and is labeled b v,j , where v represents the city and j represents the visit order. The Ising Hamiltonian can be defined by two parts: H=H A +H B (12)

[0093] H A Constraints are imposed to ensure that each city is visited once and only once. B Models the sum of the travel distances of two neighboring visited cities, where D uvis the travel distance between city u and city v. Since each travel distance is calculated twice, the total is halved. A and B are H A and H B The coefficient of , which determines the impact of constraints and travel distance on the overall Hamiltonian. In order to balance the validity of the solution and the quality of the solution, B is set to 1 and A is set to max(D uv ). In addition, city 1 is always selected as the first city to visit and thus the number of spins required is reduced to (N-1) 2 This can be understood by fixing 2N-1 spins to represent city 1 and visit order 1. This has a lasting effect on the other spins and only modifies the local field term h, which in this embodiment is added in the digital domain and does not change the remaining (N-1) 2 The coupling strength between the spins. Problem example generation and optimal solution

[0094] All the max cut problem examples used in this article are generated by the numpy random module in the python environment with the default random seed. For the unweighted max cut problem, a predefined number of "0" and "1" are first given, and shuffled by the numpy.random.shuffle function to ensure specific density. For the weighted max cut problem, the connection strength is assigned by randomly selecting a 16-bit integer (in the range of 0 to 65,535) using the numpy.random.randint function. Since the running time of the exhaustive search exceeds 10 for a single problem of problem size N>50 5 s, and continues to scale exponentially with the problem size. The optimal solution used in this paper is to run SA for a sufficiently long time, i.e., 5Nlog(N) update cycles (one update cycle means updating all spins once and corresponds to N 2 1000 trials are performed for each problem and the best solution is selected among them to ensure a high confidence in reaching the true optimal solution.

Claims

1. A method for combinatorial optimization analysis by parallel annealing, It is characterized in that include: providing a plurality of analog memristor crossbar switches having a plurality of nonvolatile two-terminal memory elements to generate an encoded matrix-represented solution to the optimization problem; arranging the analog memristor crossbars into one or more arrays having many-to-many connectivity; Identify intermediate spin states; determining a gradient of the Ising Hamiltonian based on the spin configuration across each element of the array; as well as The intermediate spin states are updated based on the gradient of the Ising Hamiltonian across each element of the array.

2. The method according to claim 1, It is characterized in that Determining the intermediate spin state comprises: Representing the spin configuration as discrete values ​​via an Ising model; and The intermediate spin states are characterized using simulation variables.

3. The method according to claim 2, It is characterized in that The Ising coupling between each spin pair in the Ising model is stored as an analog conductance value in the memristor crossbar array.

4. The method according to claim 1, It is characterized in that The memristor crossbar is further arranged to have a configuration of transistors in the memristor crossbar array.

5. The method according to claim 4, It is characterized in that The memristor crossbar array includes a single-transistor single-memristor configuration.

6. The method according to claim 1, It is characterized in that The memristor crossbar array is further connected to: Driver; multiplexer; Transimpedance amplifier; analog-to-digital converters; and Digital logic.

Images