Program, data processing device, and data processing method

Abstract

To reduce calculation cost.SOLUTION: A storage section 11 stores values of a plurality of discrete variables included in an evaluation function and values of a plurality of local fields for a plurality of respective replicas, and storage sections 12b1 to 12bR provided corresponding to the plurality of respective replicas store the values of the plurality of corresponding discrete variables and the values of the plurality of corresponding local fields. A processing section 12 repeats processing of updating any value of the plurality of discrete variables, the values of the evaluation function, and the values of the plurality of local fields, on the basis of a temperature to be set and the values of the plurality of local fields stored in the storage sections 12b1 to 12bR for the plurality of respective replicas, performs re-sampling processing of a population annealing method every time the number of repetition reaches a predetermined number, reads the values of the plurality of discrete variables and the values of the plurality of local fields for a first replica from the storage section 11 when a second replica which is a duplicate of the first replica is generated, and stores them in a storage section 12bi provided corresponding to the second replica.SELECTED DRAWING: Figure 5

Description

[Technical Field]

[0001] The present invention relates to a program, a data processing device, and a data processing method. [Background technology]

[0002] Markov chain Monte Carlo (MCMC) algorithms are widely used in discrete optimization, machine learning, medical applications, and statistical physics. Furthermore, MCMC algorithms are used in pseudo-annealing, annealed importance sampling, parallel tempering (also known as replica exchange), and population annealing (see, for example, Non-Patent Documents 1-6). Of these, population annealing is suitable for searching for solutions to large-scale combinatorial optimization problems using parallel hardware (see, for example, Non-Patent Documents 7 and 8).

[0003] In population annealing, annealing (gradually lowering the temperature) is performed, similar to the pseudo-annealing method. However, based on the solution search results for each of the multiple replicas of the Boltzmann machine at predetermined intervals, a weight (evaluation value) is determined for each replica. Then, replicas with low evaluation values ​​are eliminated, and replicas with high evaluation values ​​are duplicated (split), and the search is repeated. This process of eliminating and duplicating replicas according to their evaluation values ​​is called resampling.

[0004] Furthermore, the MCMC algorithm can be used to sample the Boltzmann distribution of a Boltzmann machine at low temperatures for the purpose of solving optimization problems or spin glass simulations (see, for example, Non-Patent Document 9). One of the main factors that negatively impacts the performance of a Boltzmann machine is the rejection of proposed state transitions (hereinafter sometimes referred to as MCMC moves). In particular, when sampling from a rough distribution, many of the proposed MCMC moves are rejected. In this case, a considerable amount of computation time is wasted because the state does not change.

[0005] To address this problem, two parallel MCMC algorithms have been proposed: the uniform selection MCMC algorithm (see, for example, Non-Patent Document 1) and the jump MCMC algorithm (see, for example, Non-Patent Document 10). Instead of proposing a single MCMC move and accepting or rejecting it, these algorithms evaluate several candidate MCMC moves in parallel to increase the likelihood that at least one MCMC move will be accepted. When implemented in parallel hardware, these algorithms show significant speed improvements, especially at low temperatures where acceptance probabilities are low. [Prior art documents] [Non-patent literature]

[0006] [Non-Patent Document 1] M. Aramon et al., “Physics-inspired optimization for quadratic unconstrained problems using a digital annealer”, Frontiers in Physics, vol. 7, Article 48, 2019 [Non-Patent Document 2] M. Bagherbeik et al., “A permutational boltzmann machine with parallel tempering for solving combinatorial optimization problems”, In International Conference on Parallel Problem Solving from Nature, pp.317-331, Springer, 2020 [Non-Patent Document 3] “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 101-1 671-680,

Outdoor Tools 4

Direct Environment 5

Outdoor Configuration 6

Direct Environment 7

Outdoor Track 8

[0007] When using population annealing to search for solutions to large-scale combinatorial optimization problems, a problem arises in that the computational cost becomes high because a large number of replicas are used.

[0008] In one aspect, the present invention aims to provide a program, a data processing device, and a data processing method that can reduce computational costs. [Means for solving the problem]

[0009] In one embodiment, for each of the multiple replicas of the Boltzmann machine, the values ​​of multiple discrete variables included in the evaluation function of the Boltzmann machine transformed from a combinatorial optimization problem, and the values ​​of multiple local fields representing the amount of change in the value of the evaluation function when each of the values ​​of the multiple discrete variables changes, are stored in a first memory unit. The corresponding values ​​of the multiple discrete variables and the values ​​of the multiple local fields are stored in a second memory unit provided corresponding to each of the multiple replicas. For each of the multiple replicas, based on the set temperature parameter value and the values ​​of the multiple local fields stored in the second memory unit, the values ​​of any of the multiple discrete variables, the value of the evaluation function, and the multiple local fields are calculated. A program is provided that causes a computer to perform a process of repeatedly updating the values ​​of a plurality of discrete variables and the evaluation function, and each time the number of iterations of the updating process reaches a predetermined number of times, a resampling process of the population annealing method is performed based on the values ​​of the plurality of discrete variables of the plurality of replicas and the value of the evaluation function, and if a second replica, which is a copy of the first replica among the plurality of replicas, is generated by the resampling process, the program reads the values ​​of the plurality of discrete variables and the plurality of local field values ​​of the first replica from the first storage unit or the second storage unit provided in correspondence with the first replica, and stores them in the second storage unit provided in correspondence with the second replica.

[0010] In one embodiment, a data processing device is provided. In one embodiment, a data processing method is provided. [Effects of the Invention]

[0011] In one respect, the present invention can reduce the computational cost when searching for solutions to combinatorial optimization problems using population annealing. [Brief explanation of the drawing]

[0012] [Figure 1] This figure shows examples of the Boltzmann distribution at various temperatures. [Figure 2]This figure shows an example of an MCMC algorithm for Boltzmann machines. [Figure 3] This figure shows an example of the overall algorithm for the population annealing method. [Figure 4] This figure shows an example of a resampling algorithm. [Figure 5] This figure shows an example of a data processing device of this embodiment that searches for solutions to combinatorial optimization problems using the population annealing method. [Figure 6] This flowchart shows an example of the procedure for searching for a solution using the population annealing method. [Figure 7] This flowchart shows an example of the resampling procedure. [Figure 8] This figure shows an example of the acceptance probability of an MCMC simulation in the maximum cut problem. [Figure 9] This figure shows an example of a parallel shrinking tree. [Figure 10] This figure shows an example of a sequential selection MCMC algorithm. [Figure 11] This figure shows an example of a parallel minimum shrinking tree. [Figure 12] This flowchart shows an example of the procedure for a sequential selection MCMC algorithm. [Figure 13] This figure shows an example implementation using a GPU. [Figure 14] This figure shows an example of thread processing within a thread block. [Figure 15] This figure shows examples of test results for sample energy histograms when uniform selection MCMC algorithms, sequential selection MCMC algorithms, and jump MCMC algorithms are applied. [Figure 16] This figure shows examples of the average error calculation results when applying the uniform selection MCMC algorithm, the jump MCMC algorithm, and the sequential selection MCMC algorithm. [Figure 17] This figure shows an example of a cluster. [Figure 18]This figure shows an example of a clustering algorithm. [Figure 19] This is a flowchart illustrating an example of the clustering procedure. [Figure 20] This figure shows an example of the weight coefficient matrix before and after clustering. [Figure 21] This figure shows an example of a clustered sequential selection MCMC algorithm. [Figure 22] This flowchart shows an example of the procedure for a clustered sequential selection MCMC algorithm. [Figure 23] This flowchart shows an example of the procedure for parallel update processing within a cluster. [Figure 24] This figure shows the simulation results of the maximum cut problem using clustered and non-clustered sequential selection MCMC algorithms. [Figure 25] This figure shows the simulation results for the performance of SA, PT, and PA in QAP. [Figure 26] This figure shows the relationship between temperature ladder size and success probability in PT (Progressive Testing). [Figure 27] This figure shows the simulation results for the performance of SA, PT, and PA in the maximum cut problem (G33). [Figure 28] This figure shows the simulation results for the performance of SA, PT, and PA in the maximum cut problem (G53). [Figure 29] This figure shows the simulation results for the performance of SA, PT, and PA in TSP. [Figure 30] This figure shows an example of the PBBM algorithm. [Figure 31] This figure shows the relationship between execution time and population size when performing 105 MCMC iterations per replica. [Figure 32] This figure shows the peak number of MCMC iterations that can be executed per second across the entire population. [Figure 33]This figure shows the calculated speed improvement compared to other instances of PBBM (GPU implementation). [Figure 34] This figure shows the results of a comparison of the accuracy of Type 1 solvers. [Figure 35] This figure shows the results of a comparison of the accuracy of Type 2 solvers. [Figure 36] This figure shows an example of computer hardware, which is an example of a data processing device. [Modes for carrying out the invention]

[0013] The embodiments for carrying out the invention will be described below with reference to the drawings. The data processing device of this embodiment searches for a solution by converting a combinatorial optimization problem into an evaluation function of a Boltzmann machine. As the evaluation function, an evaluation function (also called an energy function) of QUBO (Quadratic Unconstrained Binary Optimization) as expressed by equation (1) can be used.

[0014]

number

[0015] In equation (1), the bolded x is a vector of N discrete variables (binary variables with values ​​of 0 or 1). This vector is also called the state vector and represents the state of the QUBO model. Therefore, in the following, x may simply be referred to as the state. ij x is a discrete variable. i and x j w is a coefficient that represents the weight or strength of the bond between them (hereinafter referred to as the weight coefficient). ij If = 0, x i and x j This indicates that there is no connection between N b i These are the bias coefficients corresponding to each discrete variable.

[0016] Searching for a solution to a combinatorial optimization problem is equivalent to searching for the state of the QUBO model where the value of the evaluation function, as shown in equation (1), is minimized. The value of the evaluation function corresponds to energy. By changing the sign of the evaluation function, it is also possible to search for the state where the value of the evaluation function is maximized.

[0017] The QUBO model is equivalent to the Ising model and can be easily converted between them (see, for example, Andrew Lucas, “Ising formulations of many NP problems”, Frontiers in Physics, vol.2, article 5, Feb 2014).

[0018] The QUBO model can be used to formulate various real-world problems. Furthermore, the QUBO model is well-suited to hardware and is ideal for accelerating solution search using hardware.

[0019] Furthermore, the data processing device of this embodiment implements the functionality of a Boltzmann machine. The Boltzmann machine will be described below. (Boltzmann machine) A Boltzmann machine is a fully connected recurrent neural network that can use the energy function of QUBO. Each Boltzmann machine has a Boltzmann distribution that represents the probability of each state, which is represented by the state vector described above. The Boltzmann distribution can be expressed by equation (2) below.

[0020]

number

[0021] In equation (2), β is 1 / T (where T represents temperature), and is called the inverse temperature of the Boltzmann machine, Z β The partition function is represented by the following equation (3).

[0022]

number

[0023] The values of the partition function at different temperatures are used in spin glass simulations in statistical physics. However, since the calculation of Z β has a complexity on the order of 2 N , it is difficult to calculate when dealing with thousands of discrete variables. Therefore, a method of using the Monte Carlo algorithm to extract samples from the Boltzmann distribution and estimate the partition function is considered.

[0024] Sampling from the Boltzmann distribution is also a technique used to search for solutions to combinatorial optimization problems. As the temperature decreases (β increases), the state of the global minimum in the Boltzmann distribution begins to dominate. Therefore, the state of the global minimum has a higher probability than any other state and is more likely to be observed in samples obtained from the Boltzmann distribution.

[0025] Figure 1 is a diagram showing examples of Boltzmann distributions at various temperatures. In Figure 1, as the energy function of Equation (3), E(x)=e -0.3x sin4x + e -0.2x sin2x is used. The horizontal axis represents the value of x, and the vertical axis represents the Boltzmann probability.

[0026] The distribution becomes flatter at higher temperatures, and the maxima and minima are expanded at lower temperatures. At the lowest temperature, the state of the global minimum (the global maximum of the Boltzmann distribution) dominates the Boltzmann distribution.

[0027] When sampling from the Boltzmann distribution, in the MCMC algorithm, the value of a discrete variable can be inverted (hereinafter referred to as flipped), and it is possible to move from one state to another state one Hamming distance away. This flip is accepted with a probability expressed by the following Equation (4) according to the acceptance criterion of the Metropolis method.

[0028]

Equation

[0029] In equation (4), ΔE j (x) is ΔE j (x) = E(x') - E(x), which represents the energy increment due to the transition of the above states, and can be expressed by the following equation (5).

[0030]

number

[0031] The calculation of equation (5) has a complexity of the order of N, 2 It is less complex than directly calculating the energy, with a complexity of the order of magnitude. Figure 2 shows an example of an MCMC algorithm for a Boltzmann machine.

[0032] The inputs are state x, energy E (initial value), inverse temperature β, weight coefficient matrix w, and bias coefficient matrix b. The output is state x and energy E. First, the discrete variable x of identification number j (corresponding to a neuron in a recurrent neural network) j A flip of the value of is proposed, and ΔE is given by equation (5). j The following is calculated. The proposal of the discrete variables to be flipped is done randomly or sequentially. Then, the probability p in equation (4) j This is calculated, and a uniform random number r between 0 and 1 is generated. Here, p j >In the case of r, x j is 1-x j It is updated to E+ΔE j It will be updated.

[0033] Furthermore, if the MCMC algorithm satisfies the equilibrium conditions shown in equation (6) below, the equilibrium distribution π β It is guaranteed that it will converge to (x).

[0034]

number

[0035] In equation (6), p(x|x') represents the probability of proposing a move to x', accepting this move, and moving from x to x'. This is a sufficient condition to guarantee convergence. However, many actual MCMC algorithms satisfy the detailed equilibrium condition as shown in equation (7) below. The detailed equilibrium condition is a special case of the equilibrium condition.

[0036]

number

[0037] Detailed equilibrium conditions are not necessary, but they guarantee convergence. In the QUBO model, a simple method to satisfy the detailed balance condition is to randomly propose discrete variables whose values ​​are flipped. This method will be referred to as the random method below. In this case, the detailed balance condition can be expressed as shown in equation (8) below.

[0038]

number

[0039] In equation (8), ψ(x) is the set of all N states that are 1 Hamming distance away from x. In the QUBO model, one approach is to propose a method (sequential method) that flips the values ​​of discrete variables sequentially rather than randomly. In other words, in the k-th iteration, a flip of the value of the r-th discrete variable is proposed, where k ≡ r-1 (mod N).

[0040] Although this method violates the detailed equilibrium condition, convergence is guaranteed because the equilibrium condition is satisfied (see, for example, R. Ren and G. Orkoulas, “Acceleration of markov chain monte carlo simulations through sequential updating,” The Journal of Chemical Physics, vol. 124, no. 6, p. 064109, 2006). Furthermore, for example, A. Barzegar et al., “Optimization of population annealing monte carlo for large-scale spin-glass simulations,” Phys. Rev. E, vol. 98, p. 053308, Nov 2018, shows that the sequential method is more efficient than the random method, especially as the problem size increases. Moreover, the sequential method is easier to implement and computationally faster than the random method, which uses a random number generator, because it does not require a random number generator at the proposal stage.

[0041] (Population annealing method) In this embodiment of the data processing device, the population annealing method is used to search for solutions to combinatorial optimization problems. The population annealing method will be described below.

[0042] Population annealing method uses high temperature (β start At an initial inverse temperature of approximately 0, sampling by MCMC is initiated using a set of R replicas of the Boltzmann machine, based on the Boltzmann distribution. After a predetermined number of sampling cycles (θ cycles), a resampling process is performed, and the inverse temperature is set to Δβ = (β end (end inverse temperature)-β start ) / (N T It increases by -1) N Trepresents the number of temperature parameters (β in this example). During the resampling stage, high-energy replicas, including those trapped in local minima, are removed, and low-energy replicas are duplicated instead. The weight (evaluation) representing the importance of the i-th replica is defined based on the energy of that replica, as shown in equation (9) below.

[0043]

number

[0044] The evaluation value of equation (9) is normalized as shown in equation (10) below.

[0045]

number

[0046] And the evaluation value τ of equation (10) (i) Using this, the i-th replica is given by the following equation (11): (i) It is duplicated multiple times.

[0047]

number

[0048] After changing β to β'=β+Δβ, the population size (number of replicas) fluctuates slightly according to the following equation (12).

[0049]

number

[0050] Figure 3 shows an example of the overall algorithm for the population annealing method. The inputs are the weight coefficient matrix w, the bias coefficient matrix b, the replica number R, and the starting inverse temperature β. start , ending inverse temperature β end, a predetermined number of resampling cycles θ, and the number of temperature parameters N T The output is the state of each replica and the energy of each replica.

[0051] First, β to β start (Row 1) is set, and Δβ is (β end -β start ) / (N T -1) is set (row 2), R β R is set as (line 3). Then, state x is initialized randomly (line 4), and the initial value of E is calculated by equation (1) (line 5).

[0052] Then, from t=1 to N T The process from line 7 to line 14 is repeated until i=1. Also, from R β The processing of lines 8 to 12 is performed in parallel until u=1. Furthermore, the MCMC algorithm shown in Figure 2 is executed from u=1 to θ.

[0053] When u=θ is reached, the evaluation value W is set for each replica as described above. (i) The calculation is performed, and then the resampling process is carried out. Figure 4 shows an example of a resampling algorithm.

[0054] First, a vector variable x represents the initial state of the replica newly generated by the replication process. new , a vector variable E to represent the initial energy of the replica newly generated by the replication process. new It is initialized (line 17). Then, from i=1, R β The processing of lines 19 to 20 is performed in parallel until the following: Lines 19 to 20 perform processing based on equations (10) and (11), and rand()>τ (i) -c (i) If that is the case, c (i) is c (i) It is updated to +1 (line 20).

[0055] After that, j is initialized to 0 (line 22), then from i=1 to R βThe processing from lines 24 to 27 is performed until t is 1 to c (i) The process from lines 25 to 26 is repeated until the line 25 is created. Line 25 represents the initial state of x, which is the j-th replica newly generated by the replication. (j) new The state of the i-th replica x (i) This is substituted. Also, E represents the initial energy of the j-th replica. (j) new The energy of the i-th replica E (i) This is substituted. In line 26, j is incremented by 1.

[0056] i is R β When it reaches state x, new E is substituted into the energy E. new The value is substituted, and j is substituted for β. Then β is updated to β + Δβ. Figure 5 shows an example of a data processing device according to this embodiment that searches for a solution to a combinatorial optimization problem using the population annealing method.

[0057] The data processing device 10 is, for example, a computer and has a storage unit 11 and a processing unit 12. The memory unit 11 is, for example, a volatile memory device such as a DRAM (Dynamic Random Access Memory), or a non-volatile memory device such as an HDD (Hard Disk Drive) or flash memory. The memory unit 11 may also be HBM (High Bandwidth Memory). Furthermore, the memory unit 11 may include electronic circuits such as SRAM (Static Random Access Memory) registers.

[0058] Memory unit 11 calculates x1 to x for each of the multiple replicas of the Boltzmann machine. N Store the value of x and the values ​​of multiple local fields. N x1~x are, for example, the N discrete variables included in the evaluation function of the Boltzmann machine shown in equation (1). Multiple local fields are x1~x NThis represents the change in the value of the evaluation function when each of the values ​​changes (flips). i Local field h i (1≦i≦N) can be expressed by the following equation (13).

[0059]

number

[0060] Furthermore, the energy increment shown in equation (5) can be expressed using the local field in the following equation (14).

[0061]

number

[0062] Note that equation (5) is x j The value of x represents the energy increment when flipped, but equation (14) is x i The value of represents the energy increment when flipped. In equation (14), 2x i -1 is x i If = 0, then -1, x i When = 1, it is +1, therefore ΔE i (x) is the local field h i to, x i The value is obtained by multiplying the value by the sign corresponding to the value of .

[0063] Figure 5 shows x1~x for the i-th replica out of R replicas. N is x (i) h1~h N ga h (i) It is expressed as follows: The memory unit 11 contains the weight coefficient matrix w, which is information about the evaluation function, and the energy E for each replica. (1) ~E (R) The memory unit 11 also stores the bias coefficient matrix b, the replica number R, and the inverse temperature β. start , ending inverse temperature β end , a predetermined number of resampling cycles θ, and the number of temperature parameters N TIt may be remembered.

[0064] The processing unit 12 can be implemented by a hardware processor such as a CPU (Central Processing Unit), GPU (Graphics Processing Unit), or DSP (Digital Signal Processor). Alternatively, the processing unit 12 may be implemented by an electronic circuit such as an ASIC (Application Specific Integrated Circuit) or FPGA. The processing unit 12, for example, executes a program stored in the memory unit 11 to cause the data processing unit 10 to search for a solution. The processing unit 12 may also be a collection of multiple processors.

[0065] The processing unit 12 has multiple replica processing units, each processing a replica of the Boltzmann machine. Figure 5 shows R replica processing units 12a1 to 12aR that process R replicas. The replica processing units 12a1 to 12aR can be implemented, for example, by R processors or R processor cores. Note that there may be more than R replica processing units 12a1 to 12aR to account for the possibility of increasing the number of replicas due to resampling.

[0066] Furthermore, each of the replica processing units 12a1 to 12aR is provided with a memory unit. For example, the replica processing unit 12a1 is provided with a memory unit 12b1, the replica processing unit 12ai is provided with a memory unit 12bi, and the replica processing unit 12aR is provided with a memory unit 12bR.

[0067] Memory units 12b1 to 12bR are used that have a smaller storage capacity than memory unit 11 but can be accessed at a faster speed than memory unit 11. For example, volatile memory devices such as SRAM registers (or cache memory) are used. Note that each of memory units 12b1 to 12bR may be a storage area in a single memory device.

[0068] Each of the memory units 12b1 to 12bR stores the values ​​of multiple corresponding discrete variables and the values ​​of multiple local fields. For example, when the replica processing unit 12a1 processes the first replica, the memory unit 12b1 stores x (1) h (1) When x is stored and the replica processing unit 12aR processes the R-th replica, the storage unit 12bR contains x (R) h (R) This is stored in memory.

[0069] Since the weight coefficient matrix w is accessed in each iteration of the MCMC update process to update the local field, the memory units 12b1 to 12bR may store the weight coefficient matrix w. However, the weight coefficient matrix w is usually too large to fit into an SRAM register or cache memory, so in the example in Figure 5, it is stored in memory unit 11.

[0070] The memory units 12b1 to 12bR may be located outside the processing unit 12. The replica processing units 12a1 to 12aR determine the value of the temperature parameter (hereinafter simply referred to as temperature) set for each of the multiple replicas, and h stored in the storage units 12b1 to 12bR. (1) ~h (R) Based on the values, the update process using MCMC is repeated. The update process is performed based on an MCMC algorithm, for example, as shown in Figure 2, and updates the value of one of the discrete variables and the energy E that accompanies it. Furthermore, updates to multiple local fields are also performed in conjunction with the update of one of the discrete variables. Each time a discrete variable's value is updated, the corresponding discrete variable and local field values ​​stored in the memory unit 11 may be updated, or they may be updated all at once when performing the resampling process.

[0071] The processing unit 12 performs a resampling process based on the values ​​of multiple discrete variables and the energy E of multiple replicas each time the number of iterations of the update process reaches a predetermined number. The resampling process is performed based on an algorithm such as the one shown in Figure 4. Also, when a second replica, which is a copy of the first replica among the multiple replicas, is generated, the processing unit 12 reads the values ​​of multiple discrete variables and multiple local field values ​​of the first replica from the storage unit 11 and stores (copies) them in the storage unit provided for the replica processing unit that processes the second replica among the replica processing units 12a1 to 12aR.

[0072] For example, consider the case where, during the resampling process, the i-th replica is deleted and a copy of the 1st replica is generated. In this case, for example, as shown in Figure 5, x (1) h (1) This is read out and stored in the memory unit 12bi of the replica processing unit 12a1 that was processing the i-th replica.

[0073] In the above case, the processing unit 12 may read the values ​​of multiple discrete variables and multiple local field values ​​of the first replica from the storage unit 12b1 instead of the storage unit 11, and store them in the storage unit provided for the replica processing unit that processes the second replica.

[0074] By performing this process, it becomes unnecessary to recalculate the local field for each replicated model after each resampling operation. This reduces computational cost, increases computational efficiency, and improves solution performance, even when searching for solutions to large-scale combinatorial optimization problems involving a large number of replicas.

[0075] Figure 6 is a flowchart illustrating an example of the solution search process using the population annealing method. Figure 6 shows the flow of the search process based on the overall algorithm of the population annealing method shown in Figure 3.

[0076] First, various parameters are set (step S10). In the process of step S10, the weight coefficient matrix w, the bias coefficient matrix b, the number of replicas R, the initial inverse temperature β start , the final inverse temperature β end , a predetermined number of times θ indicating the resampling period, the number of temperature parameters N T etc. are set.

[0077] After that, the processing unit 12 determines the initial operation parameters (step S11). For example, the processing unit 12 determines, as the initial operation parameters, β = β start , Δβ=(β end - β start ) / (N T - 1), R β = R.

[0078] Next, the processing unit 12 sets the initial state (step S12). The processing unit 12 determines the initial state, for example, by randomly determining the values of x1~x N . Further, in the process of step S12, the processing unit 12 calculates the initial local field corresponding to the determined initial state from Equation (13). Also, the processing unit 12 calculates the initial energy corresponding to the determined initial state from Equation (1). The initial state and the initial local field are stored in the storage units 11, 12b1~12bR as the initial values of the states x (1) ~x (R) and the local fields h (1) ~h (R) . The initial energy is stored in the storage unit 11 as the initial values of the energies E (1) ~E (R) of each replica.

[0079] After that, the temperature loop (steps S13~S21) is performed while incrementing t by one from t = 1 to t < N T . Also, the replica loop (steps S14~S19) is performed while incrementing i by one from i = 1 to i < R β . Further, the iteration loop (steps S15~S17) is performed while incrementing u by one from u = 1 to u < θ.

[0080] In step S16, the processing unit 12 performs an update process using MCMC based on an MCMC algorithm, for example, as shown in Figure 2. When u=θ is reached, the processing unit 12 calculates the evaluation value W represented by the above equation (9). (i) Calculate (Step S18). i=R s If this is reached, the processing unit 12 performs a resampling process (step S20). t=N T The search process terminates when this condition is reached.

[0081] The data processing device 10 may output the search results (for example, the minimum energy value and the state x when the minimum energy value was obtained) when the search is completed. The search results may be displayed on a display device connected to the data processing device 10, for example, or transmitted to an external device.

[0082] Although steps S14 to S19 in Figure 6 show an example where processing is performed sequentially for each replica, the replica processing units 12a1 to 12aR shown in Figure 5 can execute steps S15 to S18 in parallel for R replicas.

[0083] Figure 7 is a flowchart showing an example of the resampling process. First, the replica loop (steps S30-S34) starts from i=1 to i <R s During this time, i is incremented by one each time.

[0084] In the process of step S31, the processing unit 12 evaluates the evaluation value τ represented by the above formula (10). (i) And the number of copies c, as expressed by the aforementioned equation (11). (i) Calculate. Subsequently, processing unit 12 performs rand()>τ (i) -c (i) Determine whether or not (step S32), rand()>τ (i) -c (i) If it is determined that this is the case, c (i) to c (i) Update to +1 (Step S33).

[0085] Processing unit 12 performs rand()>τ (i) -c (i) If it is determined that it is not, or if after processing in step S33, i is R β If the value has not been reached, increment i by 1 and return to step S31.

[0086] i is R β When this is reached, the processing unit 12 generates a vector variable x to represent the initial state, initial energy, and initial local field of the replica newly generated by the replication. new ,E new ,h new Initialize the variable and initialize the variable j to 0 (step S35).

[0087] Subsequently, the replica loop (steps S36-S41) runs from i=1 to i <R s During this time, i is incremented by one each time. Also, the replication loop (steps S37-S40) runs from t=1 to t <c (i) During this time, t is increased by one each time.

[0088] In step S38, the processing unit 12 performs a replica copy. In step S38, the processing unit 12 performs a replica copy that represents the initial state of the j-th replica newly generated by the copy. (j) new The state x of the i-th replica read from memory unit 11 (i) Substitute this value. Also, the processing unit 12 sets E, which represents the initial energy of the j-th replica. (j) new Then, the energy E of the i-th replica read from memory unit 11 (i) Substitute the value. Furthermore, the processing unit 12 generates the initial local field of the j-th replica h (j) new The local field h of the i-th replica read from memory unit 11 (i) Substitute this value.

[0089] Subsequently, the processing unit 12 substitutes j+1 for j (step S39), and t becomes c (i)If it has not reached the target value, increment t by 1 and return to step S38. t is c (i) If it reaches R, the processing unit 12 will β If it has not reached the target value, increment i by 1 and return to step S37.

[0090] i is R β When this is reached, the processing unit 12 sets the state x for the new one or more replicas generated by the replication. new Substitute E into the energy E. new Substitute h into the local field h new Substitute (step S42). The process in step S42 is to store the state x=x for the new replica in the memory unit 11. new Energy E=E new , local field h=h new This is stored. Also, the memory unit of the replica processing unit that processes new replicas (one of memory units 12b1 to 12bR) stores the state x=x new , local field h=h new It is stored (copied).

[0091] In steps S38 and S42, the values ​​of multiple discrete variables and multiple local field values ​​of the first replica are read from the memory unit 11, and stored in the memory unit of the replica processing unit that processes the second replica among the replica processing units 12a1 to 12aR.

[0092] Subsequently, the processing unit 12 updates β to β+Δβ (step S43) and returns to the process of step S21 in Figure 6. Although steps S30 to S34 in Figure 7 show an example where processing is performed sequentially for each replica, the replica processing units 12a1 to 12aR shown in Figure 5 can execute steps S31 to S33 in parallel for R replicas.

[0093] Furthermore, the processing order shown in Figures 6 and 7 is just one example, and the processing order may be changed as appropriate. (Parallel MCMC algorithm) Incidentally, in MCMC algorithms, many of the proposed MCMC migrations may be rejected during Monte Carlo simulations because they significantly increase energy or involve very low temperatures. In this case, the system remains in the same state for a long time, which is inefficient in solving combinatorial optimization problems.

[0094] Figure 8 shows an example of the acceptance probability of an MCMC simulation in the maximum cut problem. Figure 8 shows an example using Max-Cut G54, one instance of the maximum cut problem (for example, Y. ye, “Gset Max-Cut Library”, 2003, [Retrieved June 28, 2022], Internet).<URL: https: / / web.stanford.edu / yyye / yyye / Gset / > (See reference). The horizontal axis represents inverse temperature, and the vertical axis represents acceptance probability.

[0095] As is clear from Figure 8, the acceptance probability at inverse temperatures of 2.0 or higher (more than half of the entire simulation) is less than 10%. This means that more than 90% of the iterations in the MCMC simulation are wasted.

[0096] One technique to avoid rejection of proposed MCMC movement is the parallel MCMC algorithm. In the parallel MCMC algorithm, instead of proposing a flip of the value of one discrete variable (a single MCMC movement) and deciding whether to accept it or not, the acceptance probability of the value flip is evaluated in parallel for all discrete variables (this stage is called the trial stage). Then, one value of the discrete variable is flipped based on the acceptance probability. This method accelerates convergence to a target distribution in sampling problems or the detection of a global minimum state in combinatorial optimization problems.

[0097] Two examples of parallel MCMC algorithms are as follows: (Uniform Selection MCMC Algorithm) After the trial stage, the selection stage is performed. In the selection stage, the uniform selection MCMC algorithm generates random numbers in parallel for each discrete variable, and it is determined whether a value flip is potentially acceptable according to equation (4). If there are multiple discrete variables for which a value flip is acceptable, one of those discrete variables is randomly selected with a certain probability, and the value of the selected discrete variable is flipped.

[0098] This method does not guarantee convergence to an equilibrium distribution (see, for example, Non-Patent Document 9). Furthermore, in practice, randomly proposing a value flip for a single discrete variable and deciding whether it is acceptable is different from evaluating value flips for all discrete variables in parallel and randomly selecting one of the discrete variables whose value flip is accepted.

[0099] To implement the above selection stage in hardware that performs parallel processing, a parallel shrinking tree can be used, which compares pairs of discrete variables for which a value flip has been accepted and selects one discrete variable (see, for example, Non-Patent Document 8). When using a parallel shrinking tree, one discrete variable is randomly selected from the pair of discrete variables for which a value flip has been accepted.

[0100] Figure 9 shows an example of a parallel shrinking tree. Figure 9 shows an example of a parallel shrinking tree in a Boltzmann machine containing six discrete variables (x1 to x6). The discrete variables whose values ​​have been flipped are indicated by black circles. In the example in Figure 9, x1, x2, or x4 is selected, but it is not randomly selected from among x1, x2, and x4. The probabilities of x1, x2, and x4 being selected are 0.25, 0.25, and 0.5, respectively.

[0101] Therefore, implementations using parallel shrinking trees may not achieve fair and uniform selection, potentially resulting in an erroneous distribution. (Jump MCMC algorithm) The jump MCMC algorithm was developed to solve the convergence problem of the uniform selection MCMC algorithm described above (see, for example, Non-Patent Document 9). After the trial stage, the i-th discrete variable x i The probability of flipping the value of is calculated using equation (4) above. The probability for each of the N discrete variables is normalized as shown in equation (15) below.

[0102]

number

[0103] And the j-th discrete variable is p j The value is selected to be flipped with a certain probability. This reduces the rejection rate of MCMC move proposals in conventional MCMC algorithms. Since p(x|x)=0, x j Before the state x changes to x' due to a flip in the value of , the multiplicity M(x) of state x is calculated. M(x) indicates the number of iterations the state would remain in if the jump MCMC algorithm were not used. The multiplicity is a random variable from the probability distribution shown in equation (16) below.

[0104]

number

[0105] In equation (16), α(x) represents the probability of escaping from state x, and is given by equation (17) below.

[0106]

number

[0107] As a result, instead of normal samples, the jump MCMC algorithm generates weighted samples, where the weight of each sample is equal to its multiplicity. The expected value of the multiplicity of state x can be calculated as shown in equation (18) below.

[0108]

number

[0109] While the jump MCMC algorithm solves the convergence problem in the uniform selection MCMC algorithm, calculating the multiplicity of each state leads to additional overhead. Furthermore, the jump MCMC algorithm uses a random method that randomly proposes discrete variables to flip values, making it less efficient than the sequential method.

[0110] The sequential method has been applied to sparsely connected spin glass models and has been shown to be faster than the random method (see Non-Patent Documents 10 and 11). However, the sequential method has not been applied as a parallel MCMC algorithm to fully connected Boltzmann machines.

[0111] The data processing device 10 of this embodiment can, for example, apply the following sequential selection MCMC algorithm. Figure 10 shows an example of a sequential selection MCMC algorithm.

[0112] The inputs are the state x, energy E, inverse temperature β, weight coefficient matrix w, bias coefficient matrix b, and the identification number k of the neuron (discrete variable) whose value was last updated. The outputs are the state x and energy E.

[0113] From i=1 to N, the processing from row 2 to row 12 is performed in parallel. In row 2, ΔE is calculated based on equation (5) above. i This is calculated, and in row 3, based on the above formula (4), p i The following is calculated. Note that if the local field represented by equation (13) is used, ΔE is calculated based on the aforementioned equation (14). i This is calculated.

[0114] In line 4, a uniform random number between 0 and 1 is generated as r, and in line 5, an integer flag value for the discrete variable with identification number i is generated, f iIt is initialized to =2N+1, where N is the number of discrete variables in the Boltzmann machine (or Ising model).

[0115] After that, p i If >r, then the processing in lines 7 to 11 is performed. In the processing in lines 7 to 11, if i > k, then f i =i is set, otherwise f i It is set to =i+N. In the processing of line 14, the parallel minimum shrinking tree described later is used based on the obtained flag value to determine the discrete variable to be updated (represented by identification number j).

[0116] After that, f i If <2N+1, then x j is 1-x j It was updated to E + ΔE j It is updated to this. Furthermore, if the local field represented by equation (13) is used, the local field is also updated. As described above, the sequential selection MCMC algorithm, like other parallel MCMC algorithms, first runs a trial stage in parallel for all discrete variables. After the trial stage, the identification number of the discrete variable whose value flip was accepted first (the next smallest identification number after k) is selected from among the identification numbers after k. To accelerate this selection, a parallel minimum reduction tree algorithm can be used (see, for example, M. Harris, “Optimizing Parallel Reduction in CUDA”, NVIDIA Developer Technology, 2007).

[0117] Figure 11 shows an example of a parallel minimum shrinkage tree. Figure 11 shows an example of a parallel minimum-shrinking tree in a Boltzmann machine (N=6) containing six discrete variables (x1 to x6). Discrete variables for which value flipping was accepted are indicated by black circles. Discrete variables for which value flipping was not accepted are indicated by white circles. The numbers inside the circles represent the flag value f mentioned above. i That is the case.

[0118] In the selection of discrete variables to be flipped using a parallel minimum reduction tree, as in the case of using a parallel reduction tree, pairs of discrete variables are compared and one discrete variable is selected. However, when using a parallel minimum reduction tree, among the pairs of discrete variables, the one with a smaller flag value f i is selected.

[0119] The flag value f assigned to the discrete variable with identification number i i can be expressed by the following equation (19).

[0120]

Equation

[0121] As a result, for a discrete variable whose identification number i is greater than k and whose value flip is accepted, a flag value f with a small value is assigned, and for a discrete variable whose identification number i is less than or equal to k and whose value flip is accepted, a flag value f with a large value to which N is added to i as a penalty is assigned. Therefore, it is guaranteed that the flag value f with the minimum value belongs to the discrete variable whose value flip was first accepted after the previously selected discrete variable. The discrete variable for which f i = 2N + 1 is not selected. i is assigned. i is assigned. Therefore, it is guaranteed that the flag value f with the minimum value belongs to the discrete variable whose value flip was first accepted after the previously selected discrete variable. For f i = 2N + 1, the discrete variable is not selected.

[0122] As shown in Fig. 11, initially k = 0. If the value flips of x1, x2, and x4 are accepted in the first trial stage, the flag values f of x1, x2, and x4 i are such that f i = i. The flag values f of x3, x5, and x6 for which the value flip was not accepted i are such that f i = 2N + 1 = 13. In this case, as shown in Fig. 11, x1 is selected. As a result, k is updated to 1.

[0123] If the value flips of x4 and x6 are accepted in the next trial stage, the flag values f of x4 and x6 i are such that fi becomes = i. The flag values f of x1, x2, x3, and x5 for which the flip of the value was not accepted i is f i becomes = 2N + 1 = 13. In this case, as shown in FIG. 11, x4 is selected. As a result, k is updated to 4.

[0124] Hereinafter, the procedure of the sequential selection MCMC algorithm by the data processing apparatus 10 of the present embodiment will be summarized in a flowchart. In the data processing apparatus 10 that executes population annealing, each of the replica processing units 12a1 to 12aR of the processing unit 12 can perform the following processing in parallel.

[0125] FIG. 12 is a flowchart showing an example of the procedure of the sequential selection MCMC algorithm. In FIG. 12, the flow of processing based on the sequential selection MCMC algorithm shown in FIG. 10 is shown.

[0126] First, the replica processing units 12a1 to 12aR perform a trial loop (steps S50 to S58) while incrementing i one by one from i = 1 to i < N. In the trial loop, the replica processing units 12a1 to 12aR calculate ΔE i based on the above formula (14) (step S51), and calculate p i based on the above formula (4) (step S52). Thereafter, the replica processing units 12a1 to 12aR initialize the flag value f i to f i = 2N + 1 (step S53), and determine whether p i > r (a uniform random number from 0 to 1) (step S54).

[0127] p i If it is determined that > r, the replica processing units 12a1 to 12aR determine whether i > k (step S55). If the processing unit 12 determines that i > k, f i is set to = i (step S56), and if it is determined that i > k is not true, f iSet =i+N (step S57). Note that if the aforementioned population annealing is performed, the flag value f set for each replica is set. i Since it is used in each trial, it is stored in memory units 12b1 to 12bR, which can be accessed at a faster speed than memory unit 11, in Figure 5.

[0128] In the process of step S54, p i If it is determined that it is not r, or after processing in steps S56 and S57, if i has not reached N, the replica processing units 12a1 to 12aR increment i by 1 and return to the processing in step S51.

[0129] When i reaches N, the replica processing units 12a1 to 12aR use the aforementioned parallel minimum shrinking tree to update x j Determine (step S59). When using a parallel minimum shrinking tree, j is the minimum flag value f i This will be the identification number.

[0130] Subsequently, the replica processing units 12a1~12aR perform f i Determine whether or not < 2N+1 (step S60). i If it is determined that <2N+1, the replica processing units 12a1 to 12aR perform the update process (step S61). In the update process, x j is 1-x j It was updated to E + ΔE j It will be updated to x. j As the value is updated, the local field represented by equation (13) is also updated. This completes one iteration of the sequential selection MCMC.

[0131] In Figure 12, the processes in steps S50 to S58 may be executed in parallel by each of the replica processing units 12a1 to 12aR for multiple i values. Furthermore, the processing order shown in Figure 12 is just one example, and the processing order may be changed as appropriate.

[0132] (Example test) The following shows the results of tests conducted to compare the performance of applying a sequential selection MCMC algorithm, a uniform selection MCMC algorithm, and a jump MCMC algorithm when performing population annealing.

[0133] In the following test example, a GPU with 80 streaming multiprocessors (SMs) was used as the processing unit 12. Each SM has 64 processor cores capable of executing instructions in parallel. The replica processing units 12a1 to 12aR, as shown in Figure 5, can be implemented using multiple thread blocks. Each thread block contains multiple threads. The compiler automatically assigns each thread block to an SM and each thread to a processor core.

[0134] Figure 13 shows an example of implementation using a GPU. In Figure 13, elements that are the same as those shown in Figure 5 are denoted by the same reference numerals. Figure 13 shows an example where thread blocks 1 to R correspond to replica processing units 12a1 to 12aR. In this implementation example, the memory units 12b1 to 12bR shown in Figure 5 are on-chip shared memory accessible to all threads within a thread block. Communication between threads within a thread block is possible via this shared memory.

[0135] When assigning the processing of a group of ν discrete variables to each thread, the number of threads is N / ν. If the maximum number of threads in a single thread block is 1024, and N > 1024, then ν > 1. In the following experimental example, ν = 16 is used.

[0136] In each iteration of parallel MCMC processing, each thread has a probability p that it will flip one of the discrete variables assigned to it. j Evaluate the flag f j This generates ΔE based on equation (5). j (x) is calculated, and its ΔE jUsing (x), by equation (4), p j The value is calculated. Flag f j This is stored in the shared memory mentioned above (for example, in thread block 1, memory unit 12b1). The thread also executes a parallel shrinking tree algorithm. When applying the sequential selection MCMC algorithm, the parallel shrinking tree is a parallel minimum shrinking tree as shown in Figure 11; when applying the uniform selection MCMC algorithm, a tree as shown in Figure 9 is used. When applying the jump MCMC algorithm, probability p j Instead, normalized probabilities, as shown in equation (15), are used.

[0137] Figure 14 shows an example of thread processing within a thread block. Figure 14 shows an example of thread processing in the replica processing unit 12ai (corresponding to thread block i) which processes the i-th replica. In the parallel trial stage, N / ν threads execute trials in parallel, and flags f1~f N It is stored in shared memory (storage unit 12bi).

[0138] The flag is either 0 or 1 for the uniform selection MCMC algorithm, a probability as shown in equation (15) for the jump MCMC algorithm, and calculated from equation (19) for the sequential selection MCMC algorithm.

[0139] Subsequently, the shrinking tree stage is executed. The vertically aligned threads in the diagram are synchronized. In each thread, one of two flags is selected. First, the convergence of the three MCMC algorithms described above was tested using a large-scale problem in which the Boltzmann distribution is precisely known. Next, instances of the maximum cut problem were computed using each of these MCMC algorithms.

[0140] The two-dimensional Ising model with bimodal disorder can be solved exactly and thus serves as an appropriate benchmark for testing each of the above MCMC algorithms. In this problem, N discrete variables with values of -1 or +1 and N×N weight coefficients with values of -1, 0, or +1 are used. When the N variables are arranged in a two-dimensional grid, the weight coefficients between each neuron and its four adjacent neurons (left, right, above, and below) are randomly set to -1 or +1, and all other weight coefficients are set to 0.

[0141] Instances of the two-dimensional Ising model with 1024 spins are generated according to a known energy distribution. The energy distribution is represented by the following equation (20).

[0142]

Equation

[0143] Here, an instance has two ground states: a state where all values of the 1024 discrete variables are +1 and a state where all values are -1. The above three MCMC algorithms were run for N T θ = 1.024×10 5 iterations using replicas with T = 1 / β = 10 and R = 640. Then, 30 samples were obtained from each replica, resulting in a total of m ≈ 19200 samples.

[0144] Figure 15 shows an example of the test results of the energy histograms of samples when the uniform selection MCMC algorithm, the sequential selection MCMC algorithm, and the jump MCMC algorithm are applied. The horizontal axis represents energy, and the vertical axis represents probability.

[0145] In Figure 15, the exact distribution is shown together with the energy histograms when each MCMC algorithm is applied. As shown in Figure 15, the energy histogram of the sample obtained when the uniform selection MCMC algorithm was applied deviated significantly from the accurate distribution. In contrast, the energy histogram of the sample obtained when the sequential selection MCMC algorithm used by the data processing device 10 of this embodiment was applied closely matched the accurate distribution, similar to the case when the jump MCMC algorithm was applied. Therefore, it can be confirmed that the sampling was performed correctly.

[0146] Next, we present the test results of calculating 60 instances of the maximum cut problem obtained from the aforementioned "Gset Max-Cut Library" using the three MCMC algorithms described above. The maximum cut problem involves coloring the nodes of a graph G with black and white in such a way that the sum of the weights of the edges connecting the black and white nodes is maximized. The energy function of QUBO for the maximum cut problem can be expressed by the following equation (21).

[0147]

number

[0148] In equation (21), e represents a set of edges in graph G, and e ij This represents the weight of the edge connecting nodes i and j. Furthermore, equation (21) can be converted to the form of equation (1).

[0149] To compare the performance of the three aforementioned MCMC algorithms for the maximum cut problem, the best value (best cut value) and the average value (average cut value) were used. For each instance, each MCMC algorithm was run 100 times. The parameters shown in Table 1 below were used.

[0150] [Table 1]

[0151] The simulation results are shown in Tables 2 to 4.

[0152] [Table 2]

[0153] [Table 3]

[0154] [Table 4]

[0155] In Tables 2 to 4, G1 to G60 are instance names ("Name"). "Sequential" represents the simulation results when using the sequential selection MCMC algorithm, "Jump" represents the simulation results when using the jump MCMC algorithm, and "Uniform" represents the simulation results when using the uniform selection MCMC algorithm. The value in parentheses in the best cut value results is the number of times that best cut value was obtained (hit count).

[0156] The Best Known Solution (BKS) can be obtained, for example, from the following references 1-4. (Reference 1) Fuda Ma and Jin-Kao Hao, “A multiple search operator heuristic for the max-k-cut problem”, Annals of Operations Research, 248(1), pp. 365-403, 2017 (Reference 2) VP Shylo, OV Shylo, and VA Roschyn, “Solving weighted max-cut problem by global equilibrium search”, Cybernetics and Systems Analysis, Vol.48, No. 4, pp. 563-567, July, 2012 (Reference 3) Fuda Ma, Jin-Kao Hao, and Yang Wang, “An effective iterated tabu search for the maximum bisection problem”, Computers and Operations Research, 81, pp. 78-89, 2017 (Reference 4) Qinghua Wu, Yang Wang, and Zhipeng Lu, “A tabu search based hybrid evolutionary algorithm for the max-cut problem”, Applied Soft Computing, 34, pp. 827-837, 2015 The error of each MCMC algorithm in each instance is defined as shown in equation (22) below.

[0157]

number

[0158] f BKS This represents the BKS in the instance, and f BKS This represents the average cut-off value. Figure 16 shows examples of the average error calculation results when applying the uniform selection MCMC algorithm, the jump MCMC algorithm, and the sequential selection MCMC algorithm. The horizontal axis represents the error (average error) (in percent). The vertical axis represents the three MCMC algorithms.

[0159] Figure 16 shows the average error of each MCMC algorithm in 60 instances, as shown in Tables 2 to 4. As shown in Figure 16, the sequential selection algorithm has a smaller average error than the uniform selection MCMC algorithm and the jump MCMC algorithm, and is more accurate when solving the maximum cut problem.

[0160] Furthermore, when performing the cluster-based parallel update process described later, using a sequential selection algorithm is easier to implement than other algorithms. (Clustering) N discrete variables (neurons) may include groups of discrete variables that are not connected to each other.

[0161] In the following, a set of n consecutive discrete variables (neuron set) that are not connected to each other is defined as C = {x a ,x a+1 ,…,x a+n This is denoted as}. The fact that n consecutive discrete variables are not connected to each other can be expressed using weight coefficients in the following equation (23).

[0162]

number

[0163] In equation (23), i and j are two discrete variables x in C. i ,x j This is the identification number. C is also called an unconnected cluster, but hereafter it will simply be referred to as a cluster. From equation (5), x i Even if you flip the value of ΔE j It is clear that this does not affect (x). Therefore, a flip in the value of one discrete variable within a cluster does not affect the probability of flipping the values ​​of other discrete variables within the cluster.

[0164] One of the main advantages of applying a sequential selection MCMC algorithm is that, in the parallel trial stage where the probability of flipping the values ​​of all discrete variables is determined, all discrete variables whose values ​​have been accepted can be updated in parallel within the same cluster in a single trial. In uniform selection MCMC and jump MCMC algorithms, the discrete variables whose values ​​are flipped are proposed randomly, so the discrete variables in consecutive iterations may not necessarily be in the same cluster, and therefore cannot be updated together in parallel.

[0165] In problems that can be represented by a QUBO model with a weight coefficient matrix w and a bias coefficient matrix b, a graph coloring problem with respect to w can be used to detect clusters. This is because discrete variables that are not connected to each other should be assigned the same color.

[0166] Figure 17 shows an example of a cluster. In the example in Figure 17, x1 to x4 are assigned the same color, x5 to x7 are assigned the same color, and x8 is assigned a different color from x1 to x7. In other words, x1 to x4 belong to the first cluster, x5 to x7 belong to the second cluster, and x8 belong to the third cluster.

[0167] The data processing device 10 of this embodiment can cluster discrete variables by solving a graph coloring problem as a preprocessing step using population annealing with a sequential selection MCMC algorithm, and can detect the above-mentioned discrete variable groups (clusters).

[0168] The graph coloring problem is defined in QUBO form using a preprocessing weight coefficient matrix (represented by w with a hat symbol) and bias coefficients (represented by b with a hat symbol), as shown in equation (24) below.

[0169]

number

[0170] The weight coefficient matrix represented by equation (24) above is the combined (w) in the energy function of the QUBO model. ij This penalizes all pairs of discrete variables (where ≠ 0). In each iteration of the MCMC algorithm, clusters of unconnected discrete variables are detected by solving a preprocessing QUBO. The weight coefficients for these clusters are then removed from the weight coefficient matrix, and the same process is repeated until all clusters are detected. After the desired proportion of discrete variables have been assigned to clusters, the remaining discrete variables are assigned to separate clusters (each discrete variable is assigned to a single cluster).

[0171] Removing the weight coefficients related to the clusters from the weight coefficient matrix in each iteration and saving the new weight coefficient matrix would incur additional overhead. Therefore, the data processing device 10 only needs to change the bias coefficients of the discrete variables belonging to the clusters to a sufficiently large number P. This ensures that these discrete variables do not exist within the clusters in the next iteration.

[0172] According to equation (5), the maximum energy increment due to a value flip can be expressed by the following equation (25).

[0173]

number

[0174] Therefore, by selecting a P such that P > 2N, we can consider the value of the discrete variable to be 0 in the next iteration. Figure 18 shows an example of a clustering algorithm.

[0175] The input to the algorithm shown in Figure 18 includes a preprocessing weight coefficient matrix (N × N real numbers) and bias coefficients (N real numbers) as shown in equation (24), as well as the proportion f (values ​​from 0 to 1) of discrete variables (neurons) to be colored. The output is the colored result (cluster classification result) c (a set of N natural numbers).

[0176] First, k is initialized to 0 (line 1), and while k < fN, the processes in lines 3 to 8 are performed. In line 3, population annealing is performed using the weight coefficient matrix and bias coefficient represented by Equation (24).

[0177] In the example of FIG. 18, as parameters used in population annealing, the number of replicas R = 640 and the starting inverse temperature β start = 0.5, the ending inverse temperature β end = 4 are used. Further, as parameters used in population annealing, a predetermined number θ indicating the period of resampling and the number of temperature parameters N T are used. N T θ is the number of iterations of the iterative process in the MCMC algorithm. N T θ is made sufficiently small so as to be negligible compared to the number of iterations when solving the original QUBO problem in order to reduce the overhead in the preprocessing. In the experimental example described later, N T = 10 and θ = N are used.

[0178] From i = 1 to N, the processes in lines 5 to 6 are performed in parallel. In line 5, b i is updated to Px i , and in line 6, k is updated to k + x i . When i reaches N, P is incremented by 1 (line 8). Depending on the value of P, it is determined to which color each discrete variable is assigned (to which cluster it is classified).

[0179] When k ≥ fN, from i = 1 to N, the processes in lines 11 to 14 are performed. In the processes of lines 11 to 14, when b i is -1, b i is updated to P, and P is incremented by 1. When i reaches N, b - P (where c, b, and P are column vectors with N elements each) is assigned to c (line 16), and c is returned to the calling function (line 17).

[0180] The following summarizes the clustering procedure by the data processing apparatus 10 of the present embodiment in a flowchart. FIG. 19 is a flowchart showing an example of the clustering procedure. In FIG. 19, the flow of processing based on the clustering algorithm shown in FIG. 18 is shown.

[0181] First, the processing unit 12 generates a graph coloring problem (step S70). As described above, the graph coloring problem is defined in the QUBO format and is generated by obtaining the weight coefficient matrix and the bias coefficient for preprocessing based on Equation (24).

[0182] Next, the processing unit 12 initializes k to 0 (step S71) and determines whether k < fN (step S72). If the processing unit 12 determines that k < fN, it performs the processing of step S73, and if it determines that k is not < fN, it performs the processing of step S78.

[0183] In the processing of step S73, the processing unit 12 executes population annealing (denoted as PA in FIG. 19) using the weight coefficient matrix and the bias coefficient represented by Equation (24). Thereafter, the processing unit 12 performs a problem update loop (steps S74 to S76) while incrementing i one by one from i = 1 to i < N.

[0184] In the problem update loop, the processing unit 12 updates b i to Px i and updates k to k + x i (step S75). After the processing of step S75, if i has not reached N, the processing unit 12 increments i by 1 and repeats the processing of step S75.

[0185] When i reaches N, the processing unit 12 increments P by 1 and repeats the processing from step S72. When the processing unit 1 determines that k is not < fN, it performs a post-processing loop (steps S78 to S81) while incrementing i one by one from i = 1 to i < N.

[0186] In the post-processing loop, the processing unit 12 performs b i Determine whether or not it is -1 (step S79). i If it is determined that is -1, the processing unit 12 will i Update to P and increment P by 1 (step S80). b i If it is determined that i is not -1, or after processing in step S80, if i has not reached N, the processing unit 12 increments i by 1 and repeats the processing from step S79.

[0187] When i reaches N, the processing unit 12 obtains the colored result c by substituting bP (where c, b, and P are column vectors with N elements each) into c (step S82), and ends the clustering.

[0188] Note that the processing order shown in Figure 19 is just one example, and the processing order may be changed as appropriate. Figure 20 shows an example of the weight coefficient matrix before and after clustering. Figure 20 shows the weight coefficient matrix for G54, one instance of the maximum cut problem with N=1000. Weight coefficients other than 0 are shown in black, and weight coefficients of 0 are shown in white. The weight coefficient matrix after clustering is rearranged so that the identification numbers (i or j) of discrete variables within the same cluster are adjacent.

[0189] In the implementation shown in Figure 13, eight discrete variables that are unconnected to all other discrete variables are added so that clustering can be performed when ν=16. Therefore, i,j can be up to 1008. Clustering is performed using the algorithm shown in Figure 18, N T The simulation was performed with =10, θ=N=1008, and f=0.99.

[0190] As a result of the clustering, in addition to the six clusters of discrete variables 346, 250, 160, 104, 58, and 40 shown in Figure 20, three more clusters (not shown) were identified through steps S72-S77 in Figure 19. Subsequently, nine discrete variables that did not belong to any cluster were obtained. These discrete variables were recognized as different clusters through steps S78-S81 in Figure 19.

[0191] The sequential selection MCMC algorithm described above can be extended to the following algorithm (called the clustered sequential selection MCMC algorithm) using the results of the clustering described above.

[0192] Figure 21 shows an example of a clustered sequential selection MCMC algorithm. The inputs are: state x, energy E, inverse temperature β, weight coefficient matrix w, bias coefficient matrix b, and m clusters C1~C m The information is (for example, the identification numbers of the discrete variables belonging to the cluster), and the identification number k of the last updated cluster. The output is the state x and energy E.

[0193] The processing from row 1 to row 14 is the same as the processing of the sequential selection MCMC algorithm shown in Figure 10. In processing line 15, x j The identification number p of the cluster to which it belongs is determined. Then, cluster C p All x included i Regarding this, the processes in lines 17 to 20 are executed in parallel.

[0194] In the processing of lines 17 to 20, if i ≥ j, then f i If <2N+1, then x i is 1-x i It was updated to E + ΔE i It is updated to this. Furthermore, if the local field represented by equation (13) is used, the local field is also updated.

[0195] The following flowchart summarizes the procedure of the clustered sequential selection MCMC algorithm by the data processing device 10 of this embodiment. Figure 22 is a flowchart illustrating an example of the steps of a clustered sequential selection MCMC algorithm. Figure 22 shows the clustering steps and the processing flow based on the clustered sequential selection MCMC algorithm shown in Figure 21.

[0196] The processing unit 12 performs clustering as shown in Figure 19 (step S90). Then, in the data processing unit 10 that performs population annealing, each of the replica processing units 12a1 to 12aR of the processing unit 12 performs the following processing.

[0197] The replica processing units 12a1 to 12aR perform parallel trials (step S91). The parallel trials are the same as the processes in steps S50 to S58 shown in Figure 12. After that, the replica processing units 12a1 to 12aR perform the process in step S92, which is the same as the process in step S59 shown in Figure 12.

[0198] The replica processing units 12a1 to 12aR flip the value x j The identification number p of the cluster to which it belongs is obtained (step S93). Then, the replica processing units 12a1 to 12aR perform parallel update processing within the cluster (step S94), and one iteration of the clustering sequential selection MCMC is completed. From the second time onward, the process from step S91 is repeated.

[0199] Figure 23 is a flowchart illustrating an example of the procedure for parallel update processing within a cluster. Figure 23 shows the processing flow based on the clustered sequential selection MCMC algorithm shown in Figure 21.

[0200] The replica processing units 12a1 to 12aR perform an intra-cluster loop (steps S100 to S103) from i=j to i <j+N p During this time, increment i by one each time. p Cluster C p x included iIt is the number of [number].

[0201] The replica processing units 12a1 to 12aR are f i Determine whether or not it is <2N+1 (step S101). i If it is determined that <2N+1, the replica processing units 12a1 to 12aR perform the update process (step S102). In the update process, x i is 1-x i It was updated to E + ΔE i It will be updated to x. i As the expression is updated, the local field represented by equation (13) is also updated.

[0202] When i reaches j+Np, the parallel update process within the cluster terminates. In Figure 23, the processing in steps S100 to S103 is executed in parallel for multiple i units by each of the replica processing units 12a1 to 12aR.

[0203] (Example test) The following shows the results of tests conducted to verify the effectiveness of the clustered sequential selection MCMC algorithm.

[0204] Figure 24 shows the simulation results of the maximum cut problem using clustered and non-clustered sequential selection MCMC algorithms. The horizontal axis represents time (seconds), and the vertical axis represents energy.

[0205] In this example, G54, one instance of the maximum cut problem, is used. The number of nodes is N=1008. The calculation conditions for population annealing are: number of replicas R=640, β start =0, end inverse temperature β end =3, Number of temperature parameters N T = 1000, and θ = 8192, which is a predetermined number of resampling cycles. The clustering and unclustering sequential selection algorithms were each executed 5 times, and Figure 24 plots the relationship between average energy and time.

[0206] The clustered sequential selection MCMC algorithm yields the same results (same minimum energy) as the non-clustered sequential selection MCMC algorithm, but the time required to obtain those results is reduced by up to three times.

[0207] Below, we present the results of testing the speed improvement of the clustered sequential selection MCMC algorithm over the non-clustered sequential selection MCMC algorithm using several more instances. Note that for clustering purposes, N T t = 10, θ = N was used in the calculations for all instances. pp t ss t css These represent the execution times for clustering (preprocessing), (non-clustered) sequential selection MCMC algorithm, and clustered sequential selection MCMC algorithm, respectively. Speedup can be calculated using the following equation (26).

[0208]

number

[0209] All instances were run under the same computational conditions used to obtain the test results shown in Figure 24.

[0210] [Table 5]

[0211] Table 5 shows the results of testing the speed improvement of the clustered sequential selection MCMC algorithm over the non-clustered sequential selection MCMC algorithm for several instances selected from the instances shown in Tables 2 to 4. pp t ss t css The unit is seconds. As shown in Table 5, clustering resulted in a speed improvement of 1.65 to 3.17 times for all instances.

[0212] (Comparison of other MCMC algorithms with population annealing) Below, we will discuss other MCMC algorithms, namely pseudo-annealing (hereinafter referred to as "SA") and parallel tempering (hereinafter referred to as "PT"). Population annealing will be referred to as "PA" below.

[0213] Unlike the Metropolis algorithm, which uses a fixed temperature, SA uses a sufficiently high temperature (β) to prevent the algorithm from getting stuck at local minima. start Sampling from the distribution begins at β. With each MCMC iteration, the temperature is lowered according to a predefined annealing schedule. This process is completed when the final temperature (β) is reached. end >β start This process is repeated until the state of the Boltzmann machine is fixed.

[0214] SA is not originally a population-based method, but to allow for a fair comparison with other population-based algorithms such as PT and PA, it typically involves R independent runs or R replicas running in parallel. It is then expected that at least one will find the global minimum state, where R is the population size (equivalent to the number of replicas in PA and PT). SA is highly likely to remain in a local minimum until the simulation ends, as it becomes even more difficult to escape that minimum as the temperature decreases if a state gets stuck in a local minimum at a certain temperature.

[0215] In PT, R replicas of the Boltzmann machine are executed in parallel. Each replica performs MCMC sampling at different fixed temperatures selected according to a predefined temperature ladder. With each iteration of the MCMC process, we propose that the replicas at adjacent temperatures exchange temperatures. The acceptance probability (SAP) of the exchange can be expressed by equation (27) below.

[0216]

number

[0217] In equation (27), i and j represent the identification numbers of two replicas with adjacent temperatures, and E (i) represents the energy of the replica with identification number i. Temperature exchange causes replicas stuck at local minimums to rise to higher temperatures, making it easier for them to escape the minimum. However, since each replica has a different temperature, the number of replicas R directly affects the temperature ladder. β start and β end If R is fixed, increasing R reduces the temperature difference between adjacent replicas, bringing SAP closer to 1. When SAP approaches 1, exchanges are accepted regardless of the energy difference between replicas, making exchange transfer inefficient. Therefore, the optimal way to increase R is to run k independent PTs in parallel and run r replicas (so that R=kr) in each PT, not so many as to degrade performance.

[0218] Unlike PT as described above, PA's multiple replicas do not depend on temperature selection. In other words, fixed β start and β end Therefore, increasing R does not affect the temperature gap. Consequently, since R can be increased arbitrarily, PA becomes a more suitable algorithm for massively parallel hardware.

[0219] (Application to three different combinatorial optimization problems) The following will cover the Quadratic Assignment Problem (QAP), the Maximum Cut Problem, and the Traveling Salesman Problem (TSP).

[0220] QAP is a problem of assigning n facilities to n locations (assigning each facility to one location) such that the sum of the products of the distance between locations and the flow between facilities is minimized. The evaluation function (energy function) of QAP can be expressed as shown in equation (28) below.

[0221]

number

[0222] In equation (28), x ik P is a Boolean variable indicating whether the i-th facility is assigned to the k-th location. P is a value sufficient to impose a penalty for a state that violates the constraint that each facility must be assigned to one location and each location to one facility.

[0223] The maximum cut problem has already been explained using equation (21), etc. The TSP problem involves visiting n cities exactly once and returning to the starting city while minimizing the total distance traveled. The evaluation function (energy function) of the TSP can be expressed as shown in equation (29) below.

[0224]

number

[0225] In equation (29), x ik This indicates whether city i is the kth city visited, and d ij is the distance between cities i and j. P is a value sufficient to impose a penalty for a state that violates the constraints that each city has been visited once and that each city occupies a specific position in the order of visited cities.

[0226] Equations (28) and (29) above can be transformed into equation (1). In this case, x in equation (1) i The number N is n 2 To become an individual. As mentioned above, the population size of SA can be the number of replicas running in parallel. A single replica performing SA is P SA (1) With probability, the solution to the problem can be found within a certain number of iterations, P SAIf (R) is defined as the success probability of R replicas, this probability can be expressed as shown in equation (30) below. This success probability is the probability that at least one replica finds the optimal solution.

[0227]

number

[0228] Since PT is replicated k times with temperature ladder size r, the probability of success is P. PT (R) can be expressed by the following equation (31).

[0229]

number

[0230] In equation (31), P PT (r) represents the probability of success when k=1. In other words, this probability of success is the probability that at least one of the r replicas finds the optimal solution.

[0231] Unlike SA and PT, the performance of PA is measured experimentally under various problems. However, when using, for example, massively parallel hardware with a large number of processor cores as the processing unit 12, it shows higher performance than SA and PT. In other words, the success probability of PA P PA (R) is expected to be as shown in equation (32) below.

[0232]

number

[0233] In a single replica, the resampling step is neutral, so in equation (32), P PA (1) = P SA (1) (QAP test) Figure 25 shows the simulation results for the performance of SA, PT, and PA in QAP. The horizontal axis represents population size, and the vertical axis represents success probability.

[0234] Figure 26 shows the relationship between the temperature ladder size and the success probability in PT. The horizontal axis represents the temperature ladder size, and the vertical axis represents the success probability. In Figure 26, the population size is assumed to be 2560.

[0235] The instance used was esc32 obtained from QAPLIB, which is the standard QAP benchmark library (see Rainer E. Burkard, Stefan E. Karisch, and Franz Rendl, “Qaplib - a quadratic assignment problem library”, Journal of Global Optimization, 10(4), pp. 391-403, 1997).

[0236] In this case, n=32, N=1024, and the optimal solution is known. All SA, PT, and PA are 5 × 10 for each replica. 5 The MCMC iteration process is performed, and β start =0, β end It was executed with a value of 1.5.

[0237] The PT results are shown with a temperature ladder size r=16. As shown in Figure 26, this has been experimentally proven to be the optimal configuration. The success probabilities for SA and PT were measured based on results from 200 runs with R=2560. Each success probability fits equations (30) and (31). For PA, the success probability was measured with various population sizes over 200 runs.

[0238] As shown in Figure 25, SA, PT, and PA each exhibit similar performance at small population sizes. However, as the population size increases, the performance of PA improves compared to SA and PT, and the success probabilities are related as shown in equation (32).

[0239] In PA, the population size at which the success rate is 99% is 1 / 3.8th of PT and 1 / 10.2th of SA. (Exam on maximum reduction problem) Figure 27 shows the simulation results for the performance of SA, PT, and PA in the maximum cut problem (G33). The horizontal axis represents population size, and the vertical axis represents the average cut value.

[0240] Figure 28 shows the simulation results for the performance of SA, PT, and PA in the maximum cut problem (G53). The horizontal axis represents the population size, and the vertical axis represents the average ITS (Iterations To Solution), i.e., the number of MCMC iterations until a known solution (BKS) is obtained.

[0241] The two instances, G33 and G53, were obtained from the aforementioned "Gset Max-Cut Library". These instances were benchmarked according to Reference 1 mentioned above. G33 had N=2000 and G53 had N=1000. The BKS for G33 and G53 were also obtained from Reference 1.

[0242] Beta in both instances start =0.2 and β end The value was set to 3.5, and the PT temperature ladder size was set to r=32. Figure 27 shows 100 runs of SA, PT, and PA (3 × 10 5 The average cutoff value detected (with MCMC iterations) is shown. As the population size increases, the average cutoff value approaches that of BKS, but PA shows higher performance because it can more efficiently avoid the state being constrained to local minimums.

[0243] Figure 28 shows that PT performs better at small population sizes (R=64), being approximately 3.3 times faster than PA. However, as the population size increases, the average ITS of PA decreases more rapidly than that of PT, and at a population size of 320, PA is faster than PT. At a population size of 1280, PA is 2.5 times faster than PT and 10.8 times faster than SA.

[0244] (TSP exam) Figure 29 shows the simulation results for the performance of SA, PT, and PA in TSP. The horizontal axis represents the average ITS in PA, and the vertical axis represents the average ITS value.

[0245] In the example in Figure 29, the average ITS of 10 randomly generated TSP instances in 32 cities (N=1024) is shown. BKS is a value validated using a known TSP solver. In this test, all instances had R=1280, optimal temperature ladder size, and common β start and β end It was executed in [location].

[0246] As shown in Figure 29, PA outperformed SA and PT in all 10 instances, being up to 4 times faster than PT and up to 9 times faster than SA. As described above, the data processing device 10 of this embodiment solves combinatorial optimization problems using PA, which has excellent performance as described above. As explained with reference to Figure 5, the data processing device 10 of this embodiment does not require the recalculation of the local field for each replicated replica each time a resampling process is performed. As a result, even in the case of a large population size (=number of replicas) R, where the performance improvement of PA is more pronounced than that of SA or PT, computational costs are reduced, computational efficiency is increased, and the solution performance is improved.

[0247] (Verification of effectiveness by data processing device 10) The following shows the results of verifying the effectiveness of the data processing device 10. The following verification uses a data processing unit 10 implemented as shown in Figure 13. A GPU with 80 SMs is used as the processing unit 12. Each SM has 64 processor cores capable of executing instructions in parallel.

[0248] The algorithm used in this verification (hereinafter referred to as the parallel PBBM (Population-Based Boltzmann Machine) algorithm) is shown below. Figure 30 shows an example of a PBBM algorithm.

[0249] The input is the same as the PA algorithm shown in Figure 3. The output is the minimum energy found and the state at that time. The processing outline of the PBBM algorithm is as follows:

[0250] β start Starting from, R β n replicas are executed in parallel at each temperature. A parallel MCMC algorithm is used for each replica to avoid the proposed MCMC move being rejected and to achieve significant speedups with each iteration. Instead of proposing to flip the value of a single discrete variable and evaluating whether to accept the value flip, the acceptance of the value flip is evaluated in parallel for all discrete variables. Then, one of the discrete variables for which the value flip is accepted is randomly selected.

[0251] Flag vector f (i) This is used for the i-th replica and indicates whether the flip of the values ​​of each discrete variable has been accepted. That is, f j (i) =1 is the x of the i-th replica j The flip of the value of f is accepted. j (i) =0 means that the flip of that value is rejected.

[0252] The parallel reduction algorithm indicated by “Parallel_Reduction” uses the flag vector f (i) It is executed using, and the identification number = k(f k (i) It returns (if = 1 exists). The parallel minimum shrinking tree described above is used as the parallel shrinking algorithm.

[0253] x j After the value of is flipped, the local fields of each discrete variable in each replica are calculated (updated) in parallel according to equation (13). After θ iterations of MCMC at each temperature, β increases by Δβ and resampling is performed on the replicas. The copy number of each replica is calculated in parallel according to equation (11). Then, the state and local fields of each replica are replaced with new values.

[0254] Figure 31 shows 10 per replica 5 This figure shows the relationship between execution time and population size during a single MCMC iteration. The horizontal axis represents population size, and the vertical axis represents execution time.

[0255] Figure 31 shows the relationship between execution time and population size for several different Boltzmann sizes (number of discrete variables) N. When N=1024 or greater, ν=16 was used as ν in Figure 13, and when N=1024 or less, ν=8 was used.

[0256] As shown in Figure 31, for N ≤ 1024, the increase in execution time is kept to less than 27% compared to when the population size is 80, up to a population size of 640.

[0257] Figure 32 shows the peak number of MCMC iterations that can be executed per second across the entire population. The horizontal axis represents the Boltzmann size (number of discrete variables), and the vertical axis represents the number of MCMC iterations per second.

[0258] In Figure 32, for comparison, the peak count reported in Non-Patent Document 8 (which implements PT using parallel trials with an FPGA) that can handle Boltzmann sizes of N ≤ 1024 is shown as "conventional". Both PA resampling and PT exchange movement have very little overhead, and the main performance criterion is the number of iterations of the MCMC iteration process that can be executed per second.

[0259] (Comparison of performance with other solvers for various instances of the maximum cut problem) In the experimental example, an instance of the maximum cut problem from the aforementioned "Gset Max-Cut Library" is used as the benchmark.

[0260] The following shows a comparison with the maximum cut solvers in the following references A-G, where the average time to reach BKS is reported. Reference A (same as Reference 1 mentioned above): Fuda Ma and Jin-Kao Hao, “A multiple search operator heuristic for the max-k-cut problem”, Annals of Operations Research, 248(1), pp. 365-403, 2017 Reference B (same as reference 4 mentioned above): Qinghua Wu, Yang Wang, and Zhipeng Lu, “A tabu search based hybrid evolutionary algorithm for the max-cut problem”, Applied Soft Computing, 34, pp. 827-837, 2015. Reference C (same as reference 3 mentioned above): Fuda Ma, Jin-Kao Hao, and Yang Wang, “An effective iterated tabu search for the maximum bisection problem”, Computers Operations Research, 81, pp. 78-89, 2017. Reference D (same as Reference 2 mentioned above): VP Shylo, OV Shylo, and VA Roschyn, “Solving weighted max-cut problem by global equilibrium search”, Cybernetics and Systems Analysis, Vol.48, No. 4, pp. 563-567, July, 2012. Reference E: A. Yavorsky et al., “Highly parallel algorithm for the ising ground state searching problem”, arXiv:1907.05124v2 [quant-ph] 16 July, 2019 Literature F: Yu Zou and Mingjie Lin, “Massively simulating adiabatic bifurcations with fpga to solve combinatorial optimization”, FPGA '20: Proceeding of the 2020 ACM / SIGDA International Symposium on Field-Programmable Gate Arrays, February, 2020 Reference G: Chase Cook et al., “Gpu-based ising computing for solving max-cut combinatorial optimization problems”, Integration, the VLSI Journal, 69, pp. 335-344, 2019 Table 6 summarizes these solvers and their respective platforms.

[0261] [Table 6]

[0262] Table 6 shows the names of the solvers from references A to G mentioned above, as well as the algorithm, type, and platform (CPU, GPU, or FPGA) of each solver, which will be used when citing them in the experimental examples described later. In this embodiment, the solver is referred to as "PBBM". Furthermore, in order to make a fairer comparison with other solvers, a CPU-based implementation of PBBM is used in addition to a GPU-based implementation.

[0263] There are two types: Type 1 and Type 2. Type 1 solvers are those that run long enough to reach the BKS and publish their results. Type 2 solvers report the average cut value obtained through short runs (mostly less than 1 second).

[0264] To evaluate the performance of each solver, the average TTS (Time To best-known Solution) was measured for each instance after running the calculation 20 times. The calculations were performed for G1-G54 instances of the "Gset Max-Cut Library" where N ≤ 3000.

[0265] In PBBM, all instances are beta start =0.2, β end The process runs with a value of 3.5, R=640 (for GPU implementations), R=100 (for CPU implementations), and a linear β annealing schedule (an annealing schedule where β increases linearly). PBBM also uses a 1-hour timeout. Most instances resolve without problems with this annealing schedule. For the few instances where this annealing schedule did not work, β was used. end R values ​​of 2.5, 4.0, 4.5, and 5.0 were tried. Additionally, for some difficult instances that could not be resolved with R=640 (for GPU implementations) or R=100 (for CPU implementations), R was increased (increased the number of replicas).

[0266] Tables 7 and 8 show the measurement results for the average TTS.

[0267] [Table 7]

[0268] [Table 8]

[0269] If the solver fails to obtain the instance's BKS in all 20 runs, the average TTS is undefined and is represented as "X". As shown in Tables 7 and 8, the PBBM (GPU implementation) executed by the data processing device 10 of this embodiment significantly outperforms the other solvers in 51 out of 54 instances.

[0270] PBBM is a highly parallelized algorithm designed for massively parallel hardware, but even in CPU implementations, it outperforms all other solvers by 17 instances.

[0271] Figure 33 shows the calculation results of the speed improvement compared to other instances of PBBM (GPU implementation). The horizontal axis represents the magnitude of the speed improvement (geometric mean of speedup in equation (26) above), and the vertical axis represents the five instances compared to PBBM (GPU implementation).

[0272] As shown in Figure 33, the PBBM (GPU implementation) achieves at least a 43-fold speed improvement compared to other instances. Note that the average TTS measurement results shown in Tables 7 and 8 depend heavily on the programming language, compiler, and coding skills used. To address this issue, a comparison was made using the average error of each solver. The average error of each solver can be defined by the following equation (33).

[0273]

number

[0274] In equation (33), f avg (G i ) is instance G i The average cut obtained over 20 runs for f is shown below. best-known (G i ) is instance G i This shows the BKS.

[0275] Furthermore, for each solver, the number of instances that could be solved with 100% confidence, i.e., f, is calculated. best-known (G i )=f avg (G i A comparison was made of the number of instances of ). These comparison results are shown in Figure 34.

[0276] Figure 34 shows a comparison of the accuracy of Type 1 solvers. The horizontal axis represents five solvers, the left vertical axis represents the average error over 54 instances, and the right vertical axis represents the number of instances in which BKS was not obtained even once in 20 runs (number of instances that could not be solved with 100% confidence).

[0277] In Figure 34, the left bar represents the average error for each solver, and the right bar represents the number of instances. PBBM (GPU implementation) solved 3 out of 54 instances with 100% confidence, outperforming all other solvers. Furthermore, the average error of PBBM (GPU implementation) was approximately one-eighth of the average error of the other best solvers, indicating that PBBM (GPU implementation) is the most accurate solver.

[0278] Since all Type 2 solvers are hardware-accelerated, PBBM (GPU implementation) was used in this experimental example. For the Type 2 solvers MARS, SB, and SA, the average time to reach BKS and the average acquired cut value are reported for 17 instances. For each instance, the number of iterations of the MCMC process was set to the number that resulted in the shortest execution time among all other solvers. Each instance was run 10 times. Also, R=640, β start =0.2, β end The calculation was performed with a value of 3.5. The results for average time and average cut value are shown in Table 9 below.

[0279] [Table 9]

[0280] As shown in Table 9, PBBM achieved better average cut-off values ​​than all other solvers in 15 out of 17 instances, even though it ran in the minimum time across all instances.

[0281] Figure 35 shows a comparison of the accuracy of Type 2 solvers. The horizontal axis represents the error (mean error), and the vertical axis represents the four Type 2 solvers. As shown in Figure 35, PBBM (GPU implementation) is the most accurate of the Type 2 solvers, with an error approximately 1 / 2.5 that of the other solvers.

[0282] This concludes the explanation of the performance comparison between PBBM and other solvers for various instances of the maximum cut problem. The processing described above performed by the data processing device 10 in this embodiment (for example, Figures 6, 7, 12, 19, 22, and 23) can be implemented by software by having the data processing device 10 execute a program.

[0283] Programs can be stored on computer-readable storage media. Examples of storage media include magnetic disks, optical disks, magneto-optical disks, and semiconductor memory. Magnetic disks include flexible disks (FDs) and hard disk drives (HDDs). Optical disks include compact discs (CDs), CD-R (recordable) / RW (rewritable), digital versatile discs (DVDs), and DVD-R / RWs. Programs may be distributed on portable storage media. In such cases, the program may be copied from the portable storage media to other storage media and executed.

[0284] Figure 36 shows an example of computer hardware, which is an example of a data processing device. The computer 20 includes a processor 21, RAM 22, HDD 23, GPU 24, input interface 25, media reader 26, and communication interface 27. The above unit is connected to a bus.

[0285] The processor 21 functions as the processing unit 12 in Figure 5. The processor 21 is a processor such as a GPU or CPU, which includes arithmetic circuits that execute program instructions and memory circuits such as cache memory. The processor 21 loads at least a portion of the programs and data stored in the HDD 23 into the RAM 22 and executes the program. The processor 21 may have multiple processor cores to execute multiple replicas or multiple threads in parallel, for example, as shown in Figure 13. The computer 20 may also have multiple processors. A collection of multiple processors (multiprocessor) may be called a "processor".

[0286] RAM22 functions as the memory unit 11 in Figure 5. RAM22 is a volatile semiconductor memory that temporarily stores programs executed by the processor 21 and data used by the processor 21 for calculations. For example, RAM22 may be made of HBM. Note that the computer 20 may be equipped with other types of memory besides RAM22, and may be equipped with multiple types of memory.

[0287] HDD23 is a non-volatile storage device that stores software programs such as the OS (Operating System), middleware, and application software, as well as data. The programs include, for example, a program that instructs the computer 20 to perform a process to search for a solution to the aforementioned combinatorial optimization problem. The computer 20 may also be equipped with other types of storage devices such as flash memory or SSDs (Solid State Drives), and may have multiple non-volatile storage devices.

[0288] The GPU 24 outputs an image (for example, an image representing the search results for a solution to a combinatorial optimization problem) to a display 24a connected to the computer 20, according to instructions from the processor 21. The display 24a can be a CRT (Cathode Ray Tube) display, a liquid crystal display (LCD), a plasma display panel (PDP), or an organic electro-luminescence (OEL) display.

[0289] The input interface 25 acquires input signals from input devices 25a connected to the computer 20 and outputs them to the processor 21. Input devices 25a can include pointing devices such as mice, touch panels, touchpads, and trackballs, as well as keyboards, remote controllers, and button switches. Furthermore, multiple types of input devices may be connected to the computer 20.

[0290] The media reader 26 is a reading device that reads programs and data recorded on the recording medium 26a. Examples of recording media 26a include magnetic disks, optical disks, magneto-optical disks (MO), and semiconductor memory. Magnetic disks include floppy disks (FD) and hard disks (HDD). Optical disks include CDs and DVDs.

[0291] The media reader 26 copies programs and data read from the recording medium 26a to other recording media such as RAM 22 or HDD 23. The read programs are executed by the processor 21, for example. The recording medium 26a may be a portable recording medium and may be used for distributing programs and data. The recording medium 26a and HDD 23 may also be referred to as computer-readable recording media.

[0292] The communication interface 27 is connected to the network 27a and communicates with other information processing devices via the network 27a. The communication interface 27 may be a wired communication interface connected by a cable to a communication device such as a switch, or it may be a wireless communication interface connected by a wireless link to a base station.

[0293] The above describes one aspect of the program, data processing device, and data processing method of the present invention based on embodiments, but these are merely examples and the invention is not limited to those described above. [Explanation of symbols]

[0294] 10 Data Processing Devices 11,12b1~12bR Storage part 12 Processing Units 12a1~12aR Replica Processing Unit

Claims

1. The values ​​of multiple discrete variables included in the evaluation function of the Boltzmann machine transformed from a combinatorial optimization problem, and the values ​​of multiple local fields representing the change in the value of the evaluation function when each of the values ​​of the multiple discrete variables changes, are stored in the first memory unit for each of the multiple replicas of the Boltzmann machine. A second storage unit, provided for each of the plurality of replicas, stores the values ​​of the corresponding plurality of discrete variables and the values ​​of the plurality of local fields. For each of the plurality of replicas, the process of updating one of the plurality of discrete variables, the value of the evaluation function, and the value of the plurality of local fields based on the set temperature parameter value and the value of the plurality of local fields stored in the second memory unit is repeated. Each time the number of iterations of the update process reaches a predetermined number, Based on the values ​​of the multiple discrete variables of the multiple replicas and the value of the evaluation function, a resampling process of the population annealing method is performed. When the resampling process generates a second replica which is a copy of the first replica among the plurality of replicas, the values ​​of the plurality of discrete variables and the plurality of local field values ​​of the first replica are read from the first storage unit or the second storage unit provided in correspondence with the first replica, and stored in the second storage unit provided in correspondence with the second replica. A program that instructs a computer to perform a process.

2. The program according to claim 1, wherein the updating process for the plurality of replicas is executed in parallel by a plurality of replica processing units.

3. Based on the weight coefficients included in the evaluation function, which represent the strength of the bonds between the plurality of discrete variables, groups of discrete variables that are not connected to each other are detected. In the aforementioned updating process, the values ​​of the first set of discrete variables included in the discrete variable group are allowed to be updated in a single iteration. The program according to claim 1, which causes the computer to perform the processing.

4. The determination of whether or not to accept the inversion of each of the aforementioned discrete variables is performed in parallel. In the update process described above, if the inversion of the values ​​of the second set of discrete variables included in the set of discrete variables is accepted, the value of the discrete variable among the second set of discrete variables that has the next smallest identification number after the identification number of the discrete variable whose value was last inverted in the previous update process described above is inverted. The program according to claim 1, which causes the computer to perform the processing.

5. A first storage unit stores, for each of the multiple replicas of the Boltzmann machine, the values ​​of multiple discrete variables included in the evaluation function of the Boltzmann machine obtained by transforming a combinatorial optimization problem, and the values ​​of multiple local fields representing the amount of change in the value of the evaluation function when each of the values ​​of the multiple discrete variables changes. A second storage unit is provided corresponding to each of the plurality of replicas and stores the values ​​of the corresponding plurality of discrete variables and the values ​​of the plurality of local fields. For each of the plurality of replicas, the processing is repeated to update the value of one of the plurality of discrete variables, the value of the evaluation function, and the value of the plurality of local fields based on the set temperature parameter value and the value of the plurality of local fields stored in the second storage unit, and each time the number of iterations of the updating process reaches a predetermined number of times, the resampling process of the population annealing method is performed based on the value of the plurality of discrete variables and the value of the evaluation function of the plurality of replicas, and if the resampling process generates a second replica which is a copy of the first replica among the plurality of replicas, the processing unit reads the value of the plurality of discrete variables and the value of the plurality of local fields of the first replica from the first storage unit or the second storage unit provided in correspondence with the first replica and stores them in the second storage unit provided in correspondence with the second replica, A data processing device having

6. Computers The values ​​of multiple discrete variables included in the evaluation function of the Boltzmann machine transformed from a combinatorial optimization problem, and the values ​​of multiple local fields representing the change in the value of the evaluation function when each of the values ​​of the multiple discrete variables changes, are stored in the first memory unit for each of the multiple replicas of the Boltzmann machine. A second storage unit, provided for each of the plurality of replicas, stores the values ​​of the corresponding plurality of discrete variables and the values ​​of the plurality of local fields. For each of the plurality of replicas, the process of updating one of the plurality of discrete variables, the value of the evaluation function, and the value of the plurality of local fields based on the set temperature parameter value and the value of the plurality of local fields stored in the second memory unit is repeated. Each time the number of iterations of the update process reaches a predetermined number, Based on the values ​​of the multiple discrete variables of the multiple replicas and the value of the evaluation function, a resampling process of the population annealing method is performed. When the resampling process generates a second replica which is a copy of the first replica among the plurality of replicas, the values ​​of the plurality of discrete variables and the plurality of local field values ​​of the first replica are read from the first storage unit or the second storage unit provided in correspondence with the first replica, and stored in the second storage unit provided in correspondence with the second replica. Data processing method.

Images