Program, data processing device and data processing method

By dynamically adjusting the number of parallel trial bits based on search information, the method addresses inefficiencies in existing combinatorial optimization methods, enhancing solution performance and reducing unnecessary calculations.

JP7755157B2Active Publication Date: 2025-10-16FUJITSU LTD
View PDF 4 Cites 0 Cited by

Patent Information

Application Number
JP2022026770
Authority / Receiving Office
JP · JP
Patent Type
Patents
Current Assignee / Owner
Filing Date
2022-02-24
Publication Date
2025-10-16
Estimated Expiration
2042-02-24

AI Technical Summary

Technical Problem

Existing methods for solving combinatorial optimization problems using parallel trials in data processing devices face inefficiencies due to insufficient solution performance when the number of state variables for parallel trials is reduced, leading to difficulties in selecting appropriate state variables and generating effective state transitions.

Method used

A program and data processing method that dynamically adjusts the number of parallel trial bits based on search information, allowing for adaptive selection of state variables to be updated, thereby optimizing the calculation process and improving solution performance.

Benefits of technology

The method enhances the performance of solving large-scale combinatorial optimization problems by optimizing the amount of calculation and promoting appropriate state transitions, reducing unnecessary computations and improving solution quality.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure 0007755157000007
    Figure 0007755157000007
  • Figure 0007755157000008
    Figure 0007755157000008
  • Figure 0007755157000009
    Figure 0007755157000009
Patent Text Reader

Abstract

To improve solving performance of a combinatorial optimization problem.SOLUTION: A processing part 12 of a data processing device 10 for searching a solution of a combinatorial optimization problem represented by an energy function including a plurality of state variables determines whether to receive changes in respective values of a plurality of first state variables in parallel about the plurality of first state variables selected from the plurality of state variables, performs search processing for searching a solution by performing processing for changing the value of any state variable determined to receive a change in the value while changing the plurality of first state variables to be selected, specifies the number of the plurality of first state variables to be selected on the basis of a search situation of the search processing or the search information showing a search record of another combinatorial optimization problem, and repeats the search processing.SELECTED DRAWING: Figure 1
Need to check novelty before this filing date? Find Prior Art

Description

[Technical Field]

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

[0002] Data processing devices are sometimes used to solve combinatorial optimization problems. The data processing device converts the combinatorial optimization problem into an energy function of an Ising model, which is a model that represents the behavior of spin in magnetic materials, and searches for a combination of state variable values ​​included in the energy function that minimizes the value of the energy function. The combination of state variable values ​​that minimizes the value of the energy function corresponds to the ground state or optimal solution represented by the set of state variable values. Note that hereinafter, the value of the energy function may also be referred to as energy.

[0003] Methods for obtaining approximate solutions to combinatorial optimization problems in a practical amount of time include the simulated annealing (SA) method and the replica exchange method, which are based on the Markov-Chain Monte Carlo (MCMC) method.

[0004] In order to efficiently solve combinatorial optimization problems (search for solutions), it is possible to increase the parallelism of the solution search process. For example, a data processing device has been proposed that, in one trial (one Monte Carlo step process) to determine the state variables whose values ​​are to be updated, determines in parallel for multiple state variables whether or not to allow each update based on the amount of change in energy accompanying the update.

[0005] However, even if the calculation and judgment of the amount of change in energy for multiple state variables is performed in parallel and a large number of state variables are allowed to be updated, due to the principle of minimizing the Ising-type energy function using the MCMC method, only one state variable is updated in each trial. Therefore, as the problem scale increases, there is a possibility that unnecessary calculations will increase and the amount of calculation will increase.

[0006] In order to reduce the amount of unnecessary calculations, a method has been proposed in which a combinatorial optimization problem is divided into a plurality of subproblems, and the above-mentioned trials for each subproblem are performed in parallel (hereinafter referred to as partial parallel trials). [Prior art documents] [Patent documents]

[0007] [Patent Document 1] Japanese Patent Application Publication No. 2020-46997 [Patent Document 2] Patent Publication No. 2021-33341 [Patent Document 3] Patent Publication No. 2021-131695 Summary of the Invention [Problem to be solved by the invention]

[0008] In partially parallel trials, there is a possibility that sufficient solution performance may not be achieved depending on the number of state variables for which trials are performed in parallel (hereinafter referred to as the number of parallel trial bits). For example, depending on the problem, if the number of parallel trial bits is reduced, the number of state variables that are allowed to be updated may be too small, making it difficult to select appropriate state variables to be updated in order to minimize energy, and making it difficult to generate appropriate state transitions.

[0009] In one aspect, the present invention aims to provide a program, a data processing device, and a data processing method that can improve the performance of solving combinatorial optimization problems. [Means for solving the problem]

[0010] In one embodiment, there is provided a program that causes a computer to execute a process of searching for a solution to a combinatorial optimization problem represented by an energy function including a plurality of state variables, the program causing the computer to execute a search process that searches for the solution by determining in parallel whether or not to accept a change in the value of a plurality of first state variables selected from the plurality of state variables, changing the value of any state variable that is determined to accept the change in value while changing the plurality of selected first state variables, and identifying the number of the plurality of first state variables to be selected based on search information that indicates the search status of the search process or the search record of another combinatorial optimization problem, and repeating the search process.

[0011] Also provided in one embodiment is a data processing device. Also, in one embodiment, a data processing method is provided. [Effects of the Invention]

[0012] In one aspect, the present invention can improve the performance of solving combinatorial optimization problems. [Brief explanation of the drawings]

[0013] [Figure 1] FIG. 1 is a diagram illustrating a data processing device according to a first embodiment. [Figure 2] FIG. 10 illustrates an example of hardware of a data processing device according to a second embodiment. [Figure 3] FIG. 2 is a diagram illustrating an example of functions of a data processing device. [Figure 4] FIG. 2 illustrates an example of a module processing unit. [Figure 5] FIG. 10 is a diagram illustrating an example of a function of local field update in a module processing unit. [Figure 6] FIG. 10 is a diagram illustrating a first example of replica processing according to a determined group configuration. [Figure 7] FIG. 10 is a diagram illustrating a second example of replica processing according to a determined group configuration. [Figure 8] FIG. 10 is a diagram illustrating an example of pipeline processing. [Figure 9] FIG. 10 is a diagram illustrating an example of reading out weighting coefficients. [Figure 10] 10 is a flowchart illustrating an example of a processing procedure of the data processing device. [Figure 11] 10 is a flowchart illustrating an example of a procedure for collecting and recording search information. [Figure 12] 10 is a flowchart showing a first example of a procedure for determining the number P of parallel trial bits. [Figure 13] 10 is a flowchart showing a second example of the procedure for determining the number P of parallel trial bits. [Figure 14] 10 is a flowchart showing a third example of the procedure for determining the number P of parallel trial bits. [Figure 15] 10 is a flowchart showing an example of a procedure for parallel processing by four groups. DETAILED DESCRIPTION OF THE INVENTION

[0014] Hereinafter, embodiments of the invention will be described with reference to the drawings. [First embodiment] A first embodiment will be described.

[0015] FIG. 1 is a diagram illustrating a data processing device according to a first embodiment. The data processing device 10 searches for a solution to a combinatorial optimization problem using the MCMC method and outputs the searched solution. For example, the data processing device 10 uses an SA method or a replica exchange method based on the MCMC method to search for a solution. The data processing device 10 has a memory unit 11 and a processing unit 12.

[0016] The storage unit 11 may be a volatile storage device such as a random access memory (RAM) or a non-volatile storage device such as a flash memory. The storage unit 11 may include an electronic circuit such as a register. The processing unit 12 may be an electronic circuit such as a central processing unit (CPU), a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA), or a graphics processing unit (GPU). The processing unit 12 may be a processor that executes a program. The term "processor" may include a collection of multiple processors (a multiprocessor).

[0017] Combinatorial optimization problems are formulated using an Ising-type energy function, and can be replaced with the problem of minimizing the value of the energy function, for example. The energy function is sometimes called the objective function or evaluation function. The energy function includes multiple state variables. The state variables are binary variables that take on values ​​of 0 or 1. The state variables may be expressed as bits. The solution to the combinatorial optimization problem is represented by the values ​​of the multiple state variables (hereinafter sometimes referred to as a state vector). The solution that minimizes the value of the energy function represents the ground state of the Ising model and corresponds to the optimal solution to the combinatorial optimization problem. The value of the energy function is expressed as energy.

[0018] The Ising-type energy function is expressed by equation (1).

[0019]

number

[0020] The state vector x has multiple state variables as elements and represents the state of the Ising model. Equation (1) is an energy function formulated in the QUBO (Quadratic Unconstrained Binary Optimization) format. Note that for problems that involve maximizing energy, the sign of the energy function can be reversed.

[0021] The first term on the right side of equation (1) is the sum of the values ​​of the two state variables and the weighting coefficients for all combinations of two state variables that can be selected from all state variables, without omissions or duplications. The subscripts i and j are the indexes of the state variables. i is the i-th state variable. x j is the jth state variable. W ij is the weighting coefficient that indicates the strength of the connection between the i-th state variable and the j-th state variable. ij =W ji and W ii =0.

[0022] The second term on the right side of equation (1) is the sum of the products of the biases of all state variables and the values ​​of the state variables. i indicates the bias for the i-th state variable. Problem information including the weighting coefficients and biases included in the energy function is stored in the storage unit 11.

[0023] State variable x i The value of 1-x changes i Then, the state variable x i The increase in δx i =(1-x i )-x i =1-2x i Therefore, for the energy function E(x), the state variable x i The change in energy ΔE due to the change in i is expressed by equation (2).

[0024]

number

[0025] h i is called a local field and is expressed by equation (3). The local field may also be called a local field (LF).

[0026]

number

[0027] State variable x j When changes in the local field h i The change in δh i (j) is expressed by equation (4).

[0028]

number

[0029] The storage unit 11 stores local fields h corresponding to each of a plurality of state variables. i The processing unit 12 holds the state variable x j When the value of δh changes, i (j) h i By adding to the bit-inverted state, h i get.

[0030] The processing unit 12 calculates the change in energy ΔE i The state transition is as follows: i The Metropolis method or the Gibbs method is used to determine whether to accept a change in the value of ΔE. Specifically, in a neighborhood search that searches for a transition from a certain state to another state with lower energy than the current state, the processing unit 12 probabilistically accepts not only states with lower energy but also transitions to states with higher energy. For example, the probability A of accepting a change in the value of the state variable that causes ΔE is expressed by equation (5).

[0031]

number

[0032] β is the reciprocal (β = 1 / T) of T (T>0), which represents temperature, and is called the inverse temperature. The min operator indicates taking the minimum value among the arguments. The upper right side of Equation (5) corresponds to the Metropolis method. The lower right side of Equation (5) corresponds to the Gibbs method. The processing unit 12 compares a uniform random number u with A where 0 < u < 1 for a certain index i, and if u < A, it accepts the change in the value of the state variable x i and changes the value of the state variable x i If u ≥ A, the processing unit 12 does not accept the change in the value of the state variable x i and does not change the value of the state variable x i According to Equation (5), the larger the value of ΔE, the smaller A becomes. Also, the smaller β is, that is, the larger T is, the more likely a state transition with a large ΔE is to be allowed. For example, when the Metropolis method is used, the processing unit 12 may perform a transition determination using Equation (6) obtained by transforming Equation (5).

[0033]

Equation

[0034] That is, when ΔE satisfies Equation (6) for a uniform random number u (0 < u ≤ 1), the processing unit 12 accepts the change in the value of the corresponding state variable. When ΔE does not satisfy Equation (6) for a uniform random number u, the processing unit 12 does not accept the change in the value of the corresponding state variable.

[0035] In the data processing apparatus 10 of the first embodiment, the processing unit 12 determines a state variable whose value is to be changed (hereinafter referred to as the state variable to be updated) by partial parallel trials for the number of parallel trial bits. Further, the processing unit 12 has a function of changing the number of parallel trial bits.

[0036] FIG. 1 illustrates a part of the flow of the processing performed by the processing unit 12. (S1) The processing unit 12 first performs a search process, for example, based on the number of parallel trial bits P1. If the number of state variables included in the energy function is N, in the process of step S1, the processing unit 12 determines whether or not to accept a change in the value of the state variables (including the calculation process of ΔE) for x1 to x N In addition, the processing unit 12 changes the value of one of the state variables to be updated, which is one of the state variables (hereinafter referred to as update candidate state variables) that are determined to accept a change in value based on the above determination of the P1 state variables. When there are multiple update candidate state variables, one state variable is selected as the update target state variable randomly or according to a predetermined rule.

[0037] If the number of state variables to be updated is often 0, no state transition occurs, resulting in a waste of calculation time. Therefore, the processing unit 12 may always change the value of one of the P1 state variables in each partial parallel trial. Hereinafter, this method will be referred to as the rejection-free method.

[0038] When the rejection-free method is used, the processing unit 12 calculates the state variables x belonging to the P1 state variables. i Generate a uniform random number u[i] for max(0,ΔE i )+Tlog(-log(u[i])) i as the state variable to be updated. Note that the max operator indicates that the maximum value among the arguments is taken. For example, when the number of state variables to be updated is 0, the processing unit 12 may select one state variable to be updated by a rejection-free method.

[0039] The processing unit 12 searches for a solution by performing the above-mentioned process while changing the selected P1 state variables. N is divided into regions A1 to An (represented as parallel trial regions in Figure 1), each containing P1 state variables. For example, search is performed sequentially from region A1 to region An. Note that each region may contain the same state variables. Also, after searching up to region An, search may be performed again from region A1.

[0040] (S2) For example, when the above search process has been performed for a predetermined period of time, the processing unit 12 changes the number of state variables (number of parallel trial bits) selected in the partial parallel trial from P1 to P2 based on search information indicating the search status of the search process of step S1.

[0041] The search information indicating the search status may be, for example, the cumulative total of the number of update candidate state variables obtained in the search process for a predetermined period, or the cumulative total of the number of state variables whose values ​​have actually changed. N The information may be the amount of movement of the state vector (expressed as a Hamming distance) represented by the set of (a) and (b) whether or not the minimum value of the energy has been updated (or the number of updates). Note that the processing unit 12 may perform a search by specifying an appropriate number of parallel trial bits based on search information, which is a record of the number of parallel trial bits in searches performed in the past for other combinatorial optimization problems.

[0042] The above search information is stored in the storage unit 11 during the search process in step S1. In the processing of step S2, the processing unit 12 calculates the average value of the number of update candidate state variables in each partial parallel trial, for example, from the cumulative value of the number of update candidate state variables obtained in the search processing over a predetermined period. Then, for example, if the average value is smaller than a first threshold, the processing unit 12 changes P1 to P2, which is larger than P1. If there are few update candidate state variables, it is difficult to select appropriate state variables for minimizing energy as update targets, which may result in a deterioration in solution performance. Therefore, to promote appropriate state transitions and improve solution performance, the processing unit 12 increases the number of parallel trial bits as described above. If the average value is larger than a second threshold (> first threshold), the processing unit 12 changes P1 to P2, which is smaller than P1. This is because if there are too many update candidate state variables, only one state variable is selected as the update target, which increases unnecessary calculations and increases the amount of calculations.

[0043] Examples of methods for adjusting the number of parallel trial bits when using other search information will be described later (see FIGS. 13 and 14). (S3) After changing the number of parallel trial bits, the processing unit 12 performs search processing based on the number of parallel trial bits P2. The processing in step S3 is performed in the same manner as the processing in step S1 above. In the example of FIG. 1, x1 to x N shows an example where there are m (<n) regions B1 to Bm each containing P2 state variables. For example, the search is performed in order from region B1 to region Bm. Each region may contain the same state variables. Also, after the search up to region Bm, the search may be performed again from region B1.

[0044] When the search processing in step S3 is performed for a predetermined period, the processing unit 12 may perform the processing in step S2 based on the search information indicating the search status of the search processing in step S3, further change the number of parallel trial bits, and repeat the search processing.

[0045] Note that the processing unit 12 may perform the above search processing in parallel for a plurality of replicas each indicating a plurality of state variables. Examples of search processing using a plurality of replicas will be described in the second embodiment.

[0046] When the SA method is performed in the processing of steps S1 and S3, the processing unit 12 decreases the value of T, which is a parameter representing temperature, according to a predetermined temperature parameter change schedule, for example, each time the partial parallel trial is repeated a predetermined number of times. Then, the processing unit 12 outputs, for example, the state vector obtained when the partial parallel trial is repeated a predetermined number of times as the calculation result for the combinatorial optimization problem (for example, it may be displayed on a display device not shown). Note that the processing unit 12 may update the value (energy) of the energy function expressed by equation (1) each time a change occurs in the value of the state variable, and store the energy and state when the energy is the minimum so far in the storage unit 11. In this case, the processing unit 12 may output, for example, the state corresponding to the minimum energy stored after the partial parallel trial is repeated a predetermined number of times as the calculation result.

[0047] When the processing unit 12 performs the replica exchange method, the processing unit 12 performs the above steps S1 to S3 for each of multiple replicas, each of which has a different T value. While a specific example will be described later, the same number of parallel trial bits may be set for each replica, or different numbers of parallel trial bits may be set for each replica. The processing unit 12 performs replica exchange every time the partial parallel trial is repeated a predetermined number of times. For example, the processing unit 12 selects two replicas having adjacent T values ​​and exchanges the T values ​​or states between the two selected replicas with a predetermined exchange probability based on the energy difference between the replicas or the difference in T values. For example, the processing unit 12 updates the value of the energy function (energy) every time a change occurs in the value of a state variable in each replica, and stores the energy and state when the minimum energy is reached in the storage unit 11. Then, the processing unit 12 outputs, as a calculation result, the state corresponding to the minimum energy among all replicas, among the minimum energies stored in each replica after the partial parallel trial is repeated a predetermined number of times.

[0048] In the data processing device 10 of the first embodiment described above, the number of parallel trial bits (the number of state variables for which it is determined in parallel whether or not to accept a change in value) of partial parallel trials is changed based on the search information. This makes it possible to set the number of parallel trial bits according to the search situation that reflects the characteristics of the problem, optimizes the amount of calculation required to change the value of one state variable, and improves the performance of solving large-scale problems.

[0049] In addition to optimizing the amount of calculation, by changing the number of parallel trial bits as described above, when the value of a certain state variable changes, it is possible to adjust the period until the state variable becomes an acceptable state for the next change in value (becomes an update candidate). This makes it possible to avoid a situation where, when a change in the value of the state variable causes the state to escape from a local optimum, the value of the state variable changes again, causing the state to be constrained to a local optimum again.

[0050] [Second embodiment] Next, a second embodiment will be described. FIG. 2 illustrates an example of hardware of a data processing device according to the second embodiment.

[0051] The data processing device 20 is a computer that searches for a solution to a combinatorial optimization problem using the MCMC method and outputs the searched solution. The data processing device 20 has a CPU 21, a RAM 22, an HDD (Hard Disk Drive) 23, a GPU 24, an input interface 25, a media reader 26, a NIC (Network Interface Card) 27, and an accelerator card 28.

[0052] The CPU 21 is a processor that executes program instructions. The CPU 21 loads at least a portion of the program and data stored in the HDD 23 into the RAM 22 and executes the program. The CPU 21 may include multiple processor cores. The data processing device 20 may also have multiple processors. The processing described below may be executed in parallel using multiple processors or processor cores. A set of multiple processors may also be called a "multiprocessor" or simply a "processor."

[0053] The RAM 22 is a volatile semiconductor memory that temporarily stores programs executed by the CPU 21 and data used in calculations by the CPU 21. Note that the data processing device 20 may include a type of memory other than RAM, or may include multiple memories.

[0054] The HDD 23 is a nonvolatile storage device that stores software programs such as an OS (Operating System), middleware, and application software, as well as data. Note that the data processing device 20 may also include other types of storage devices, such as a flash memory or an SSD (Solid State Drive), or may include multiple nonvolatile storage devices.

[0055] The GPU 24 outputs an image to a display 101 connected to the data processing device 20 in accordance with an instruction from the CPU 21. The display 101 may be any type of display, such as a CRT (Cathode Ray Tube) display, a liquid crystal display (LCD: Liquid Crystal Display), a plasma display, or an organic EL (OEL: Organic Electro-Luminescence) display.

[0056] The input interface 25 receives an input signal from an input device 102 connected to the data processing device 20 and outputs the signal to the CPU 21. The input device 102 may be a pointing device such as a mouse, a touch panel, a touch pad, or a trackball, a keyboard, a remote controller, or a button switch. In addition, multiple types of input devices may be connected to the data processing device 20.

[0057] The medium reader 26 is a reading device that reads programs and data recorded on the recording medium 103. For example, a magnetic disk, an optical disk, a magneto-optical disk (MO: Magneto-Optical disk), a semiconductor memory, etc. can be used as the recording medium 103. Magnetic disks include flexible disks (FD: Flexible Disks) and HDDs. Optical disks include compact discs (CDs) and digital versatile discs (DVDs).

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

[0059] The NIC 27 is an interface that is connected to the network 104 and communicates with other computers via the network 104. The NIC 27 is connected to a communication device such as a switch or a router by a cable. The NIC 27 may also be a wireless communication interface.

[0060] The accelerator card 28 is a hardware accelerator that uses the MCMC method to search for a solution to a problem expressed by the Ising-type energy function of equation (1). The accelerator card 28 can be used as a sampler that samples states that follow the Boltzmann distribution at a given temperature by performing the MCMC method at a constant temperature or the replica exchange method, which exchanges the states of the Ising model between multiple temperatures. To solve combinatorial optimization problems, the accelerator card 28 performs annealing processes such as the replica exchange method or the SA method, which gradually decreases the value of T.

[0061] The SA method is a method for efficiently finding an optimal solution by sampling states that follow a Boltzmann distribution at each value of T and lowering the T used for sampling from a high temperature to a low temperature, i.e., by increasing the inverse temperature β. Even on the low temperature side, i.e., when β is large, some state change occurs, increasing the likelihood of finding a good solution even if the value of T is quickly reduced. For example, when using the SA method, the accelerator card 28 repeatedly attempts state transitions at a constant value of T a certain number of times and then reduces the value of T.

[0062] The replica exchange method is a technique in which the MCMC algorithm is run independently using multiple values ​​of T, and the value of T (or state) is appropriately exchanged for the state obtained for each value of T. A good solution can be found efficiently by searching a narrow range of the state space using MCMC at low temperatures and a wide range of the state space using MCMC at high temperatures. For example, when using the replica exchange method, the accelerator card 28 performs parallel trials of state transitions for each of multiple values ​​of T, and after each set number of trials, exchanges the value of T for the state obtained for each value of T with a predetermined exchange probability.

[0063] The accelerator card 28 includes an FPGA 28a. The FPGA 28a implements a search function in the accelerator card 28. The search function may be implemented by other types of electronic circuits, such as a GPU or an ASIC. The FPGA 28a includes a memory 28b. The memory 28b stores data such as problem information used in the search by the FPGA 28a, solutions searched by the FPGA 28a, and search information indicating the search status. The FPGA 28a may include multiple memories including the memory 28b.

[0064] The FPGA 28a is an example of the processing unit 12 of the first embodiment. The memory 28b is an example of the storage unit 11 of the first embodiment. The accelerator card 28 may have a RAM outside the FPGA 28a, and data stored in the memory 28b may be temporarily saved to the RAM in accordance with the processing of the FPGA 28a.

[0065] A hardware accelerator that searches for a solution to an Ising-type problem, such as the accelerator card 28, is sometimes called an Ising machine or a Boltzmann machine. The accelerator card 28 performs a solution search in parallel using multiple replicas. The replicas represent multiple state variables included in the energy function. In the following description, the state variables are referred to as bits. Each bit included in the energy function is associated with an integer index and is identified by the index.

[0066] FIG. 3 is a diagram illustrating an example of functions of the data processing device. The data processing device 20 has an overall control unit 30, M (M is an integer of 2 or more) modules (which can also be called circuit units) 31a1, 31a2, ..., 31aM, a search information aggregation unit 32, and a selector 33. The overall control unit 30, modules 31a1 to 31aM, search information aggregation unit 32, and selector 33 are realized using electronic circuits of FPGA 28a and memory 28b.

[0067] The overall control unit 30 controls the modules 31a1 to 31aM, the search information aggregation unit 32, and the selector 33. The overall control unit 30 also receives the search information aggregated by the search information aggregation unit 32, and determines the number of parallel trial bits P. Then, based on the determined P, the overall control unit 30 determines a group configuration of the modules 31a1 to 31aM (described later), and supplies group configuration information indicating the determined group configuration to the selector 33.

[0068] Furthermore, the overall control unit 30 updates the state vector of each replica held in the storage unit based on the flip bit index for each group supplied from the selector 33. The flip bit index is the index of the bit to be updated (hereinafter referred to as the flip bit).

[0069] Furthermore, the overall control unit 30 may update the energy of each replica by adding ΔE corresponding to the index to the energy stored in an energy storage unit that stores the energy corresponding to the current state vector of each replica. Note that in Fig. 3, the storage unit that stores the current state vector corresponding to each replica and the energy storage unit that stores the energy corresponding to the current state vector of each replica are omitted. The storage unit and the energy storage unit may be realized, for example, by a storage area of ​​the memory 28b in the FPGA 28a, or by a register.

[0070] Moreover, overall control unit 30 supplies control information, group configuration information, and flip bit information (hereinafter referred to as flip bit information) to modules 31a1 to 31aM. Flip bit information includes, for example, a flip bit index and a flip direction of the flip bit (information indicating whether the flip has occurred from 0 to 1 or from 1 to 0).

[0071] The modules 31a1 to 31aM each include a module control unit 31b1, 31b2, ..., 31bM, and a module processing unit 31c1, 31c2, ..., 31cM. Module control units 31b1 to 31bM receive control information, group configuration information, and flip bit information from overall control unit 30, and control pipeline control within modules 31a1 to 31aM, local field update processing of each replica, and the like.

[0072] The modules 31a1 to 31aM are appropriately combined and grouped into n groups (n is an integer equal to or greater than 2) based on the number of parallel trial bits P. The n groups, each including one or more modules, perform partial parallel trials for n of the multiple replicas in each unit processing period, based on the number of parallel trial bits P. The module processing units 31c1 to 31cM also send search information indicating the search status to the search information aggregation unit 32. Examples of the module processing units 31c1 to 31cM will be described later.

[0073] The search information aggregation unit 32 aggregates the search information and sends the aggregated search information to the overall control unit 30 . The selector 33 changes the selector configuration based on the group configuration information received from the overall control unit 30. Then, when there are multiple indices of update candidate bits (hereinafter referred to as flip candidate bits) for each group, the selector 33 selects one for each group in parallel. Then, the selector 33 outputs the selected index as a flip bit index and supplies it to the overall control unit 30.

[0074] Below, a case where the number of modules is M=8 will be described, but the number is not limited to this. Fig. 4 is a diagram showing an example of a module processing unit, in which the overall control unit 30, module control units 31b1 to 31bM, and search information aggregation unit 32 shown in Fig. 3 are omitted.

[0075] In the example of FIG. 4, the data processing device 20 has module processing units 31c1 to 31c8. Module processing unit 31c1 has a memory unit 40a, h calculation units 40b1-40bK, ΔE calculation units 40c1-40cK, selector 40d, and search information acquisition unit 40e. The other module processing units 31c2-31c8 have a similar configuration. That is, module processing unit 31c2 has a memory unit 41a, h calculation units 41b1-41bK, ΔE calculation units 41c1-41cK, selector 41d, and search information acquisition unit 41e. Module processing unit 31c3 has a memory unit 42a, h calculation units 42b1-42bK, ΔE calculation units 42c1-42cK, selector 42d, and search information acquisition unit 42e. Module processing unit 31c4 has a memory unit 43a, h calculation units 43b1-43bK, ΔE calculation units 43c1-43cK, selector 43d, and search information acquisition unit 43e. Module processing unit 31c5 has a memory unit 44a, h calculation units 44b1-44bK, ΔE calculation units 44c1-44cK, selector 44d, and search information acquisition unit 44e. Module processing unit 31c8 has a memory unit 47a, h calculation units 47b1-47bK, ΔE calculation units 47c1-47cK, selector 47d, and search information acquisition unit 47e. K is the number of bits handled by each of module processing units 31c1-31c8.

[0076] For example, memory units 40a-47a are realized by a plurality of memories including memory 28b in FPGA 28a. h calculation units 40b1-47bK, ΔE calculation units 40c1-47cK, selectors 40d-47d, and search information acquisition units 40e-47e are realized by electronic circuits in FPGA 28a.

[0077] In Fig. 4, the names of the h calculation units 40b1-47bK are written with the subscript n, such as "hn" calculation unit, to make it easier to understand that they correspond to the nth bit. Similarly, in Fig. 4, the names of the ΔE calculation units 44c1-44cK are written with the subscript n, such as "ΔEn" calculation unit, to make it easier to understand that they correspond to the nth bit.

[0078] For example, h calculation unit 40b1 and ΔE calculation unit 40c1 perform calculations on the first bit of the N bits, and h calculation unit 40bK and ΔE calculation unit 40cK perform calculations on the i-th bit.

[0079] As described above, the modules 31a1 to 31aM are appropriately combined and grouped based on the number of parallel trial bits P, and each group performs partial parallel trials for a certain replica.

[0080] 4 shows a case where module 31a1 is classified into group A, module 31a2 is classified into group B, modules 31a3 and 31a4 are classified into group C, and modules 31a5 to 31a8 are classified into group D. In this case, partial parallel trials with the number of parallel trial bits P=K are performed in groups A and B, partial parallel trials with the number of parallel trial bits P=K×2 are performed in group C, and partial parallel trials with the number of parallel trial bits P=K×4 are performed in group D.

[0081] The data processing device 20 performs partial parallel trials on multiple replicas in parallel using n processes (pipelines) for n groups, thereby enabling efficient use of the computational resources of the FPGA 28a. For example, in the example of FIG. 4, the data processing device 20 processes multiple replicas in parallel using four pipelines for groups A to D. In this example, the number of replicas is 16. The 16 replicas are denoted as replicas R0, R1, ..., R15.

[0082] Here, the information stored in the memory units 40a to 47a will be described. Each of the memory units 40a to 47a stores a weighting coefficient W={W γ,δ When the number of bits in the state vector is N, the total number of weight coefficients is N 2 W γ,δ =W δ,γ W γ,γ= 0. Since the processing of each replica is processing the same problem, the total number of weight coefficients stored does not change even if the number of replicas increases.

[0083] In the example of FIG. 4, the memory unit 40a stores the weighting coefficient W 1,1 ~W 1,N ,…,W i,1 ~W i,N For example, the weighting factor W 1,1 ~W 1,N is used in an operation corresponding to the first bit of N bits. The total number of weighting coefficients stored in memory unit 40a is i × N. Note that when the number of bits handled by each of module processing units 31c1 to 31c8 is K, i = K.

[0084] The memory unit 41a stores the weighting coefficient W i+1,1 ~W i+1,N ,…,W j,1 ~W j,N The memory unit 42a stores the weighting coefficient W j+1,1 ~W j+1,N ,…,W k,1 ~W k,N The memory unit 43a stores the weighting coefficient W k+1,1 ~W k+1,N ,…,W l,1 ~W l,N The memory unit 44a stores the weighting coefficient W l+1,1 ~W l+1,N ,…,W m,1 ~W m,N The memory unit 47a stores the weighting coefficient W o+1,1 ~W o+1,N ,…,W N,1 ~W N,N Remember.

[0085] For example, the memory units 40a to 47a are supplied with the index of the bit whose value has changed from the module control units 31b1 to 31bM. Then, the weighting coefficient corresponding to the index is read from the memory units 40a to 47a and supplied to the h calculation units 40b1 to 47bK.

[0086] 4, when the number of groups is four, a maximum of four indexes are simultaneously supplied to the memory units 40a to 47a. As a result, a maximum of four weighting factors are simultaneously supplied to each of the h calculation units 40b1 to 47bK. The four weighting factors correspond to four replicas.

[0087] The following description will be mainly focused on h calculation unit 40b1 and ΔE calculation unit 40c1, which correspond to the first bit. The other h calculation units and ΔE calculation units have the same functions. The h calculation unit 40b1 uses the weighting coefficients read from the memory unit 40a to calculate the local field h1 for each of the four replicas processed in parallel based on equations (3) and (4). For example, the h calculation unit 40b1 has a register that holds the local field h1 previously calculated for the corresponding replica, and updates the h1 of the corresponding replica stored in the register by adding the δh1 of the corresponding replica to the corresponding h1. Note that a signal indicating the inversion direction of the bit indicated by the index to be inverted for each replica is supplied from the module control unit 31b1 to the h calculation unit 40b1. The initial value of h1 is calculated in advance using equation (3) according to b1 corresponding to the problem, and is preset in the register of the h calculation unit 40b1.

[0088] ΔE calculation unit 40c1 calculates ΔE1, which is the amount of energy change corresponding to the inversion of the own bit in one replica, based on equation (2), using the local field h1 of the replica to be processed next, which is held in h calculation unit 40b1. ΔE calculation unit 40c1 can determine the inversion direction of the own bit, for example, from the current value of the own bit of the replica. For example, if the current value of the own bit is 0, the inversion direction is from 0 to 1, and if the current value of the own bit is 1, the inversion direction is from 1 to 0. ΔE calculation unit 40c1 supplies the calculated ΔE1 to selector 40d.

[0089] The selector 40d performs a determination using equation (6) for each ΔE simultaneously supplied from the ΔE calculation units 40c1 to 40cK, and determines whether or not to invert the corresponding bit. For example, the selector 40d determines whether or not to allow inversion of the bit with index=1 for the energy change ΔE1 calculated by the ΔE calculation unit 33a1, based on equation (6). Specifically, the selector 40d determines whether or not to invert the corresponding bit for the corresponding replica, based on a comparison between -ΔE1 and thermal noise corresponding to T. The thermal noise corresponds to the product of the natural logarithm of the uniform random number u in equation (6) and T.

[0090] Furthermore, the selector 40d randomly selects one of the flip candidate bits based on a random number in accordance with equation (6) and supplies an index corresponding to the selected bit to the selector 33. Note that if there is no bit determined to be flippable, the selector 40d does not need to output an index, but when the above-mentioned rejection-free method is used, the selector 40d always outputs an index for one bit.

[0091] The selectors 41d to 47d also function in the same manner as the selector 40d with respect to the bits processed by their own modules. The search information acquisition unit 40e acquires search information in the module 31a1. For example, the search information acquisition unit 40e acquires the number of indexes (corresponding to the number of flip candidate bits) output from the selector 40d as the search information. The search information acquisition unit 40e may acquire information such as the number of flip bits in each replica and the amount of change in energy as the search information.

[0092] The search information acquisition units 41e to 47e also have the same functions as the search information acquisition unit 40e. As described above, when module 31a1 is classified into group A, module 31a2 into group B, modules 31a3 and 31a4 into group C, and modules 31a5 to 31a8 into group D, selector 33 functions as follows.

[0093] Group A consists of one module 31a1, so selector 33 includes a function (indicated as "1-1 Select") to output the index output by module processing unit 31c1 of module 31a1. Group B consists of one module 31a2, so selector 33 includes a function (indicated as "1-1 Select") to output the index output by module processing unit 31c2 of module 31a2. Group C consists of two modules 31a3 and 31a4. Therefore, selector 33 includes a function (indicated as "2-1 Select") to select and output the index output by either module processing unit 31c3 or 31c4 of modules 31a3 and 31a4. Group D consists of four modules 31a5 to 31a8. Therefore, selector 33 includes a function (indicated as "4-1 Select") to select and output the index output by either module processing unit 31c5 to 31c8 of modules 31a5 to 31a8.

[0094] When selecting one of the multiple indexes, the selector 33 randomly selects one based on a random number. The selector 33 may also be configured to preferentially select one of the indexes based on selection weight information supplied from the selectors 40d to 47d, for example. The selection weight information may be, for example, the number of flip candidate bits. In this case, the index output from the module processing unit with the largest number of flip candidate bits is preferentially selected. In addition, when the selectors 40d to 47d use the rejection-free method, the max(0, ΔE i ) + Tlog(-log(u[i])) can also be used as the selection weight information. In that case, max(0, ΔE i The index output from the module processing unit that gives the smallest value of Tlog(-log(u[i]))+Tlog(-log(u[i])) is preferentially selected.

[0095] Such a selector 33 can be realized, for example, using four 8-input, 1-output gate circuits with enable. A gate circuit that realizes "1-1 Select" has 8 inputs, of which 1 is enabled by an enable signal (for example, included in group configuration information supplied from the overall control unit 30). A gate circuit that realizes "2-1 Select" has 8 inputs, of which 2 are enabled by enable signals. A gate circuit that realizes "4-1 Select" has 8 inputs, of which 4 are enabled by enable signals. Then, the above selection process is performed.

[0096] Fig. 5 is a diagram showing an example of the function of local field update in a module processing unit. Fig. 5 shows an example of the function of local field update in module processing unit 31c1 of module 31a1. The function of local field update in the other module processing units 31c2 to 31cM is the same as the function of local field update in module processing unit 31c1.

[0097] The memory unit 40a has eight memories 40p1, 40p2, ..., 40p8 corresponding to the number of modules M=8. The memory 40p1 stores the weighting coefficient W 1,1 ~W 1,i ,W 2,1 ~W 2,i ,…,W i,1 ~W i,i The memory 40p2 stores W 1,i+1 ~W 1,j ,W 2,i+1 ~W 2,j ,…,W i,i+1 ~W i,j The memory 40p8 stores W 1,k+1 ~W 1,N ,W 2,k+1 ~W 2,N ,…,W i,k+1 ~W i,N Remember.

[0098] Each of the h calculation units 40b1 to 40bK uses up to four weighting coefficients to update local fields corresponding to up to four replica bits in parallel based on equations (3) and (4). For example, the h calculation unit 40b1 has an h holding unit r1, selectors s10, s11, s12, s13, and adders c1, c2, c3, and c4. The other h calculation units have the same functions as the h calculation unit 40b1. For example, the h calculation unit 40bK has an h holding unit ri, selectors si0, si1, si2, si3, and adders c5, c6, c7, and c8. The h calculation unit 40b1 will be described below.

[0099] The h storage unit r1 stores the local field of the local bit corresponding to each of the 16 replicas. The h storage unit r1 may be configured with a flip-flop or may be configured with four RAMs that read one word per read. The local bit in the h calculation unit 40b1 is the bit with index=1.

[0100] The selector s10 selects four of the eight weighting factors read from memories 40p1 to 40p8 and supplies each of the four selected weighting factors to one of adders c1, c2, c3, and c4. The selector s10 can be realized, for example, using four eight-input, one-output gate circuits with enable. In such a gate circuit, one of the eight inputs is enabled by an enable signal supplied from the module control unit 31b1 based on the group configuration information, and the weighting factor of the enabled input is output.

[0101] The selector s11 reads the local fields of the replicas to be updated in each group from the h-storage unit r1 and supplies them to the adders c1, c2, c3, and c4. The maximum number of local fields that the selector s11 can simultaneously read from the h-storage unit r1 is four.

[0102] Adders c1, c2, c3, and c4 add the weighting coefficients output by selector s10 to the local fields related to the replicas being processed in the four groups supplied from selector s11, thereby updating the local fields and supplying them to selector s12. As described above, the signs of the weighting coefficients can be determined, for example, by a signal indicating the bit inversion direction supplied from module control unit 31b1.

[0103] The selector s12 stores the local field of the relevant replica updated by the adders c1 to c4 in the h-holding unit r1. The selector s13 reads out the local field of its own bit in the replica to be processed next in the group A to which the module 31a1 belongs from the h holding unit r1, and supplies it to the ΔE calculation unit 40c1.

[0104] In this way, the h calculation unit 40b1 can simultaneously update the local field corresponding to index=1 for a maximum of four replicas using the selectors s10, s11, and s12 and the adders c1, c2, c3, and c4.

[0105] With the above configuration, the data processing device 20 executes up to four pipelines in parallel for 16 replicas. Next, a description will be given of an example of replica processing according to the group configuration determined by the overall control unit 30. In the following example, it is assumed that the number of stages in one pipeline, that is, the number of stages, is four.

[0106] The first stage is ΔE calculation, which is a process of calculating ΔE for each bit belonging to each group in parallel. The second stage is a flip decision, which is a process of selecting one bit to be flipped for each ΔE calculated in parallel.

[0107] The third stage is W Read, which is a process of reading out weighting coefficients from the memory units 40a to 47a. The fourth stage is the h update. The h update is a process of updating the local field for the corresponding replica based on the read weight coefficient. In parallel with the h update stage, the bit to be inverted in the corresponding replica is inverted. Therefore, the h update stage can also be called the bit update stage.

[0108] The number of pipeline stages is not limited to four. In the following description, the period during which one stage of the pipeline is processed is referred to as one step period.

[0109] 6 is a diagram showing a first example of replica processing according to the determined group configuration. In Fig. 6, M0 to M7 represent modules 31a1 to 31a8. In the following figures, modules 31a1 to 31a8 are also represented as M0 to M7.

[0110] 6, first, two of the modules 31a1 to 31a8 are combined. That is, the number of groups of the modules 31a1 to 31a8 is four. In this case, for each step period, partial parallel trials are executed in parallel for four of the replicas R0 to R15, with the number of parallel trial bits P=K×2.

[0111] The data processing device 20 starts processing a replica in a group at a timing shifted by four or more step periods in each stage so that after the h update of the replica being processed in a group, the same replica is processed in the next group. This allows ΔE calculations to be performed within each replica using local fields that reflect the previous bit update, thereby observing the principle of sequential processing of the MCMC method.

[0112] In the example of Figure 6, a replica being processed in one group is processed in the next group at a timing shifted by four step periods. Next, in the example of FIG. 6, at a certain timing, the configuration changes to one in which modules 31a1 to 31a8 are combined in groups of four. That is, the number of groups consisting of modules 31a1 to 31a8 changes from four to two. In this case, for two of replicas R0 to R15, partial parallel trials are executed in parallel with the number of parallel trial bits P = K × 4 per step period. Also, since the number of replicas processed in parallel becomes two, the data processing device 20 changes the step period, for example, from the completion of the h update of a replica being processed in a certain group until the processing of the same replica in the next group. In the example of FIG. 6, when the number of groups was four, the step period was four steps, but when the number of groups was two, the step period was changed to eight steps.

[0113] FIG. 7 illustrates a second example of replica processing according to the determined group configuration. In the example of FIG. 7, modules 31a1 to 31a8 (M0 to M7) are divided into four groups, each containing 1, 1, 2, and 4 modules. In the example of FIG. 7, partial parallel trials of each replica are performed with one of four parallel trial bit counts P. Replicas R0 to R7 are processed using one module, so the parallel trial bit count P is K. Replicas R8 to R11 are processed using two modules, so the parallel trial bit count P is K×2. Replicas R12 to R15 are processed using four modules, so the parallel trial bit count P is K×4.

[0114] The cycle in which trials are performed for all N bits varies for each replica due to differences in the number of parallel trial bits P. One trial is performed for all bits in the replica processed using the fewest modules (replicas R0 to R7 in this example), so that one trial is performed for all N bits in all replicas.

[0115] 7, replicas R0 to R7 are processed using one module in one step period. If all N bits in replicas R0 to R7 are divided and processed by eight modules, the step period in which one trial for all bits is performed is 4 (number of pipeline stages) × 8 (number of modules) = 32 step periods.

[0116] In this case, the overall control unit 30 controls the modules to which each replica is assigned to process and the group configuration of each module, as shown in Figure 7, so that each replica completes at least one trial for all N bits in the 32-step period.

[0117] FIG. 8 is a diagram illustrating an example of pipeline processing. 8, the data processing device 20 starts processing a replica at a timing shifted by four step periods in each stage so that after the h update of a replica being processed in a certain group, the same replica is processed in the next group. For example, after the h update of replica R0 in the group of module 31a8 (M7), the processing of replica R0 is started at a timing shifted by four step periods in each stage so that the processing of replica R0 is started in the group of module 31a4 (M3).

[0118] This allows ΔE calculation to be performed within each replica using a local field that reflects the previous bit update, thereby maintaining the principle of sequential processing in the MCMC method. Here, the update of the local field must be reflected in all bits of the corresponding replica. Therefore, the weighting coefficients are read out simultaneously for all bits of the four replicas. As illustrated in FIG. 5, the data processing device 20 divides the memory that holds the weighting coefficients corresponding to each group into memories 40p1 to 40p8, for example. Therefore, accesses corresponding to multiple replicas do not overlap with each other to the same memory. For example, during the step periods marked with stars in FIG. 8, the weighting coefficients are read out as follows:

[0119] FIG. 9 is a diagram showing an example of reading out the weighting coefficients. Each of the memory units 40a to 47a of the modules 31a1 to 31a8 (M0 to M7) is divided into eight memories. Each of the eight memories holds a weighting coefficient between the K bits assigned to the module itself and the K bits assigned to one of the modules 31a1 to 31a8.

[0120] For example, memory unit 40a of module 31a1 stores W0(M0), W0(M1), W0(M2), W0(M3), W0(M4), W0(M5), W0(M6), and W0(M7) in eight memories (memories 40p1 to 40p8 in FIG. 5). For example, W0(M0) is a weighting factor between the K bits assigned to module 31a1. W0(M7) is a weighting factor between the K bits assigned to module 31a1 and the K bits assigned to module 31a8.

[0121] For example, memory unit 47a of module 31a8 stores W7(M0), W7(M1), W7(M2), W7(M3), W7(M4), W7(M5), W7(M6), and W7(M7) in eight separate memories. For example, W7(M0) is a weighting factor between the K bits assigned to module 31a8 and the K bits assigned to module 31a1. W7(M7) is a weighting factor between the K bits assigned to module 31a8.

[0122] In the case of the group configuration shown in FIG. 4 (groups A to D are denoted as GA to GD in FIG. 9), W0(M0) to W7(M0) are weighting coefficients used for updating the h of each module when the bits assigned to group A are inverted. W0(M1) to W7(M1) are weighting coefficients used for updating the h of each module when the bits assigned to group B are inverted. W0(M2) to W7(M2) and W0(M3) to W7(M3) are weighting coefficients used for updating the h of each module when the bits assigned to group C are inverted. W0(M4) to W7(M4), W0(M5) to W7(M5), W0(M6) to W7(M6), and W0(M7) to W7(M7) are weighting coefficients used for updating the h of each module when the bits assigned to group D are inverted.

[0123] 8, when the bit of replica R4 processed by module 31a1(M0) belonging to group A is inverted, weighting coefficients are read from the memories holding W0(M0) to W7(M0). When the bit of replica R0 processed by module 31a2(M1) belonging to group B is inverted, weighting coefficients are read from the memories holding W0(M1) to W7(M1).

[0124] Furthermore, when the bits of replica R8 processed by modules 31a3(M2) and 31a4(M3) belonging to group C are inverted, weighting factors are read from each memory holding W0(M2) to W7(M2) or W0(M3) to W7(M3). If the inverted bits are bits assigned to module 31a4, weighting factors are read from each memory holding W0(M3) to W7(M3), as shown in FIG.

[0125] Furthermore, when the bits of replica R12 processed by modules 31a5(M4) to 31a8(M7) belonging to group D are inverted, weighting coefficients are read from the memories holding weighting coefficients for any of modules 31a5(M4) to 31a8(M7). That is, weighting coefficients are read from the memories holding any of W0(M4) to W7(M4), W0(M5) to W7(M5), W0(M6) to W7(M6), or W0(M7) to W7(M7). If the inverted bit is assigned to module 31a7, weighting coefficients are read from the memories holding W0(M6) to W7(M6), as shown in FIG.

[0126] In this way, since the processing for the four replicas is processing for bits assigned to different modules, if the memory is divided by module as shown in Figure 9, memory accesses when inverting bits in the four replicas will not overlap with the same memory (same read port). This prevents memory accesses when updating h from becoming a bottleneck and lengthening the calculation time.

[0127] Note that when updating h, the data processing device 20 may determine whether the value of the weighting coefficient is 0 or not, and may read only weighting coefficients whose value is not 0 without reading from memory those whose value is 0. This reduces the number of times weighting coefficients are read from memory. Note that in this case, the number of cycles required for reading varies depending on the proportion of weighting coefficients whose value is 0 among all weighting coefficients, and the data processing device 20 may perform control to stall the pipeline when the number of cycles is longer than a predetermined threshold.

[0128] Next, a description will be given of the processing procedure of the data processing device 20. First, a description will be given of the processing procedure for one replica. FIG. 10 is a flowchart illustrating an example of a processing procedure of the data processing device.

[0129] (S20) The overall control unit 30 of the FPGA 28a performs initial settings. For example, the initial settings include setting an initial value for the number of parallel trial bits P and initializing variables for aggregating search information. In the following example, itrnum, Csum, Fsum, Dsum, Emin, and Eminupdate are used as variables for aggregating search information.

[0130] itrnum is a variable that represents the number of iterations. Csum is a variable that represents the cumulative total number of flip candidate bits. Fsum is a variable that represents the cumulative total number of flip bits. Emin is a variable that represents the minimum energy. Dsum is a variable that represents the cumulative total of the movement amount (movement distance) of the state vector, which is expressed as the Hamming distance.

[0131] In the process of step S20, itrnum is initialized to 0, Csum to 0, Fsum to 0, and Eminupdate to 0. Emin is initialized to the maximum value that the data processing device 20 can handle, for example.

[0132] In the process of step S20, for example, the overall control unit 30 may set the problem information (weighting coefficients and biases included in the energy function, etc.) supplied to the FPGA 28a under the control of the CPU 21 in the modules 31a1 to 31aM.

[0133] (S21) The overall control unit 30 determines whether it is time to change the parallel trial bit count P. For example, it is determined every predetermined period (predetermined number of iterations) that it is time to change the parallel trial bit count P. If the overall control unit 30 determines that it is time to change the number of parallel trial bits P, it proceeds to the process of step S22, and if it determines that it is not time to change the number of parallel trial bits P, it proceeds to the process of step S23.

[0134] (S22) The overall control unit 30 performs a process to determine the number of parallel trial bits P. An example of the process in step S22 will be described later. (S23) The overall control unit 30 supplies control information, group configuration information, and flip bit information to the modules 31a1 to 31aM, causing the modules 31a1 to 31aM to execute partial parallel trial loops. The overall control unit 30 also determines the group configuration of the modules 31a1 to 31aM (described later) based on the determined P, and supplies group configuration information indicating the determined group configuration to the selector 33.

[0135] (S24) A partial parallel trial is performed by any one or a combination of the modules 31a1 to 31aM, with the number of parallel trial bits being P. In the process of step S24, DE calculation and flip determination are performed in parallel for the P bits of the replica.

[0136] (S25) The selector 33 selects a flip bit. In the process of step S25, the selector 33 selects a flip bit by selecting one of the indices of the flip candidate bits obtained as a result of the flip determination. The index of the selected flip bit (flip bit index) is supplied to the overall control unit 30.

[0137] (S26) Overall control unit 30 updates the bit in the state vector of each replica held in the memory unit that corresponds to the flip bit index supplied from selector 33. Overall control unit 30 also supplies flip bit information to modules 31a1 to 31aM. Modules 31a1 to 31aM perform h update based on the flip bit information.

[0138] (S27) The search information aggregating unit 32 aggregates and records the search information. An example of the process in step S27 will be described later. (S28) Until all bits (N bits) in the replicas have been tried, modules 31a1 to 31aM shift the area in which partial parallel trials are performed and repeat the processes from step S24 to step S27 under the control of overall control unit 30. When all bits (N bits) in the replicas have been tried, overall control unit 30 proceeds to the process of step S29.

[0139] (S29) The overall control unit 30 determines whether the search has ended. The overall control unit 30 determines that the search has ended when a predetermined search end condition is met. For example, the overall control unit 30 determines that the search has ended when the number of iterations reaches a predetermined number. If it is determined that the search has ended, the FPGA 28a ends the processing. If it is determined that the search has not ended, the processing from step S21 is repeated.

[0140] When performing the SA method, the FPGA 28a decreases the value of T according to a predetermined temperature parameter change schedule, for example, every time the partial parallel trial is repeated a predetermined number of times. When performing the replica exchange method, the FPGA 28a sets different values ​​of T to each of the multiple replicas and performs replica exchange every time the partial parallel trial is repeated a predetermined number of times. For example, the FPGA 28a selects two replicas with adjacent values ​​of T and exchanges the values ​​or states of T with a predetermined exchange probability based on the energy difference or the difference in the value of T between the replicas.

[0141] When the processing is completed, the FPGA 28a outputs the state vectors corresponding to each finally obtained replica as a solution to the CPU 21. The FPGA 28a may output the energy corresponding to each replica together with the state vector to the CPU 21. The FPGA 28a may output the solution with the lowest energy among the solutions obtained by the search as the final solution to the CPU 21. The CPU 21 may control the GPU 24 to display the solution on the display 101.

[0142] Next, an example of a procedure for tallying and recording search information by the search information tallying unit 32 will be described. Fig. 11 is a flowchart showing an example of a procedure for collecting and recording search information. Note that the search information collecting unit 32 only needs to collect search information used in the process of determining the number of parallel trial bits P, but Fig. 11 shows an example in which multiple types of search information are collected.

[0143] (S40) The search information aggregation unit 32 increments itrnum by 1 (counts up). (S41) The search information aggregation unit 32 acquires search information. In this example, the search information aggregation unit 32 acquires the number of flip candidate bits C, the presence or absence of a flip F (F = 1 if present, F = 0 if absent), the current state vector Statecur, and the current energy Ecur as search information. The number of flip candidate bits C can be acquired from modules 31a1 to 31aM, and the presence or absence of a flip F can be acquired based on whether the selector 33 outputs a flip index. If the current state vector Statecur and the current energy Ecur are stored in the memory 28b, the search information aggregation unit 32 acquires Statecur and Ecur from the memory 28b.

[0144] (S42) The search information aggregation unit 32 determines whether Ecur < Emin. If the search information aggregation unit 32 determines that Ecur < Emin, it performs the process of step S43; if it determines that Ecur ≥ Emin, it performs the process of step S44.

[0145] (S43) The search information aggregation unit 32 updates Emin with Ecur and increments Eminupdate by 1 (counts up). (S44) The search information aggregation unit 32 determines whether it is the movement amount acquisition timing. For example, the search information aggregation unit 32 determines that it is the movement amount acquisition timing when itrnum has increased by a predetermined number of times from the previous movement amount acquisition timing. If the search information aggregation unit 32 determines that it is the movement amount acquisition timing, it performs the process of step S45; if it determines that it is not the movement amount acquisition timing, it performs the process of step S47.

[0146] (S45) The search information aggregation unit 32 calculates the movement amount (Hamming distance) D between the reference state vector and the current state vector Statecur. (S46) The search information aggregation unit 32 updates the reference state vector. For example, the reference state vector is updated to Statecur.

[0147] (S47) The search information aggregating unit 32 aggregates the search information. Specifically, the search information aggregating unit 32 adds C to Csum, F to Fsum, and D to Dsum, thereby updating Csum, Fsum, and Dsum.

[0148] This completes one round of processing by the search information aggregating unit 32 to tally and record the search information. The above-described collection and recording of search information may be performed for each replica, or may be performed for all replicas at once.

[0149] Next, an example of the procedure for determining the number of parallel trial bits P by the overall control unit 30 will be described. FIG. 12 is a flowchart illustrating a first example of a procedure for determining the number P of parallel trial bits.

[0150] (S50) The overall control unit 30 calculates the average value Cave of the number of flip candidate bits. The overall control unit 30 calculates Cave by dividing Csum supplied from the search information aggregating unit 32 by itrnum (the number of iterations).

[0151] (S51) The overall control unit 30 determines whether Cave>Cthu and P>Pthl. Cthu is a first threshold value of Cave. Pthl is a lower limit value of the parallel trial bit number P (for example, K (the number of bits handled by one module)). If the overall control unit 30 determines that Cave>Cthu and P>Pthl, it performs the process of step S53, and if it determines that Cave>Cthu or P>Pthl is not true, it performs the process of step S52.

[0152] (S52) The overall control unit 30 determines whether Cave < Cthl and P < Pthu. Cthl is the second threshold value of Cave, and Cthl < Cthu. Pthu is the upper limit value of the parallel trial bit number P (for example, K × M (number of modules)). When the overall control unit 30 determines that Cave < Cthl and P < Pthu, it performs the process of step S54. When it determines that Cave < Cthl is not true or P < Pthu is not true, it performs the process of step S55.

[0153] (S53) The overall control unit 30 sets P = P - Pdec in order to reduce the parallel trial bit number P. Pdec is a value that is an integer multiple of K and is predetermined. When Cave is too large, unnecessary calculations increase and the amount of computation becomes large. Therefore, in order to suppress the amount of computation, the overall control unit 30 reduces the parallel trial bit number P.

[0154] (S54) The overall control unit 30 sets P = P + Pinc in order to increase the parallel trial bit number P. Pinc is a value that is an integer multiple of K and is predetermined. Pinc may be the same value as Pdec.

[0155] When Cave is too small, it becomes difficult to select appropriate flip candidate bits in minimizing energy, and there is a possibility that the solution performance deteriorates. Therefore, in order to promote appropriate state transitions and improve the solution performance, the parallel trial bit number P is increased as described above.

[0156] (S55) The overall control unit 30 sets the determined parallel trial bit number P in modules 31a1 to 31aM. (S56) The overall control unit 30 initializes the variables for aggregating search information and ends the process of determining the parallel trial bit number P. In the process of step S56, it is initialized to itrnum = 0, Csum = 0, Fsum = 0, Dsum = 0, Eminupdate = 0.

[0157] Figure 13 is a flowchart showing a second example of the procedure for determining the parallel trial bit number P. (S60) The overall control unit 30 calculates a flip rate Frate indicating the occurrence rate of flip bits within a predetermined period. The overall control unit 30 calculates Frate by dividing Fsum supplied from the search information aggregation unit 32 by itrnum (the number of iterations).

[0158] (S61) The overall control unit 30 determines whether Frate > Fthu and P > Pthl. Fthu is the first threshold of Frate. When the overall control unit 30 determines that Frate > Fthu and P > Pthl, it performs the process of step S63. When it determines that Frate is not > Fthu or P is not > Pthl, it performs the process of step S62.

[0159] (S62) The overall control unit 30 determines whether Frate < Fthl and P < Pthu. Fthl is the second threshold of Frate, and Fthl < Fthu. When the overall control unit 30 determines that Frate < Fthl and P < Pthu, it performs the process of step S64. When it determines that Frate is not < Fthl or P is not < Pthu, it performs the process of step S65.

[0160] (S63) The overall control unit 30 sets P = P - Pdec in order to decrease the number of parallel trial bits P. When Frate is too large, state transitions occur too frequently, deteriorating the convergence of calculations, so the solution performance may deteriorate. Therefore, the overall control unit 30 decreases the number of parallel trial bits P to suppress the magnitude of Frate.

[0161] (S64) The overall control unit 30 sets P = P + Pinc in order to increase the number of parallel trial bits P. When Frate is too small, state transitions do not occur frequently enough, so the solution performance may deteriorate. Therefore, to promote state transitions and improve the solution performance, the number of parallel trial bits P is increased as described above.

[0162] The processes of steps S65 and S66 are the same as the processes of steps S55 and S56 shown in FIG. 12, and thus the description thereof is omitted. FIG. 14 is a flowchart showing a third example of the procedure for determining the number of parallel trial bits P.

[0163] (S70) The overall control unit 30 calculates the average value Dave of the movement amount D. The overall control unit 30 calculates Dave by dividing Dsum supplied from the search information aggregation unit 32 by itrnum (the number of iterations).

[0164] (S71) The overall control unit 30 determines whether Dave > Dthu, there is no update of Emin, and P > Pthl. Dthu is the first threshold value of Dave. When the overall control unit 30 determines that Dave > Dthu, there is no update of Emin, and P > Pthl, it performs the process of step S73. When the overall control unit 30 determines that Dave ≤ Dthu, or there is an update of Emin, or P ≤ Pthl, it performs the process of step S72.

[0165] Note that the presence or absence of the update of Emin can be determined by whether Eminupdate is a value of 1 or more. (S72) The overall control unit 30 determines whether Dave < Dthl, there is no update of Emin, and P < Pthu. Dthl is the second threshold value of Frate, and Dthl < Dthu. When the overall control unit 30 determines that Dave < Dthl, there is no update of Emin, and P < Pthu, it performs the process of step S74. When the overall control unit 30 determines that Dave ≥ Dthl, or there is an update of Emin, or P ≥ Pthu, it performs the process of step S75.

[0166] (S73) The overall control unit 30 sets P=P-Pdec to reduce the number of parallel trial bits P. If Emin has not been updated even though Dave is large, many unnecessary calculations may occur, which may result in a deterioration in solution performance. For this reason, the overall control unit 30 reduces the number of parallel trial bits P to prevent the occurrence of unnecessary calculations.

[0167] (S74) The overall control unit 30 sets P = P + P inc to increase the number of parallel trial bits P. If Dave is too small and Emin has not been updated, the search range is too narrow, which may be causing a deterioration in solution performance. Therefore, to widen the search range and improve solution performance, the number of parallel trial bits P is increased as described above.

[0168] The processing in steps S75 and S76 is the same as the processing in steps S55 and S56 shown in FIG. 12, and therefore a description thereof will be omitted. The process of determining the number of parallel trial bits P as described above may be performed based on the aggregation of search information for each replica, or may be performed based on the aggregation of search information for all replicas.

[0169] Furthermore, the above three types of determination processing can be combined with each other. For example, the parallel trial bit numbers P determined by the three types of determination processing are set in the modules 31a1 to 31aM.

[0170] Based on the determined value of the number of parallel trial bits P, the overall control unit 30 may adjust P so that it is the same value for each replica (see FIG. 6) or so that P is a constant ratio for each replica (see FIG. 7). This improves the efficiency of pipeline processing.

[0171] Next, the processing procedure of the data processing device 20 will be described more specifically, taking as an example a case where parallel processing is performed by four groups. Fig. 15 is a flowchart showing an example of a procedure for parallel processing by four groups. Fig. 15 includes a more specific example of the processing of steps S23 to S28 of the procedure shown in Fig. 10 for multiple replicas. The processing for determining the number of parallel trial bits P, and the processing for aggregating and recording search information are not shown.

[0172] (S80) The overall control unit 30 performs initial settings. For example, the initial settings include setting an initial value for the number of parallel trial bits P and initializing the variables used to aggregate the search information described above. Furthermore, if multiple replicas are used, the number of replicas, the number of groups (4 in the example of FIG. 15), and the replica interval between groups are set. The replica interval is set to the number of pipeline stages or a value equal to or greater than the number of stages (4 in the example of FIG. 8 described above).

[0173] Furthermore, in the processing of step S80, the groups to which each of modules 31a1 to 31aM is assigned and the replicas to be processed by each group are set first. For example, in the example shown in FIG. 7, modules 31a1 to 31a4 (M0 to M3) are assigned to one group, and replica R12 is assigned to that group. Furthermore, modules 31a5 and 31a6 (M4 and M5) are assigned to one group, and replica R8 is assigned to that group. Furthermore, module 31a7 (M6) is assigned to one group, and replica R4 is assigned to that group, and module 31a8 (M7) is assigned to one group, and replica R0 is assigned to that group.

[0174] The overall control unit 30 may assign a replica with a higher temperature (a larger value of T, which is a parameter representing the set temperature) to a group with a smaller number of modules so that the initial value of the number of parallel trial bits P is smaller.

[0175] The following four groups are denoted as G0 to G3. (S81) The modules 31a1 to 31aM perform loop processing until one trial is performed for all bits in each replica.

[0176] (S82) The overall control unit 30 assigns modules or groups to replicas for each cycle of the replica loop. In the example shown in FIG. 7, one step period corresponds to one cycle of the replica loop. The replicas assigned to each module change with each cycle of the replica loop. The modules to which the same replica is assigned change with each four cycles of the replica loop. For example, in the first cycle, replica R12 is assigned to modules 31a1 to 31a4 (M0 to M3), and replica R8 is assigned to modules 31a5 and 31a6 (M4 and M5). Also in the first cycle, replica R4 is assigned to module 31a7 (M6), and replica R0 is assigned to module 31a8 (M7). After four cycles, replica R12 is assigned to modules 31a5 to 31a8 (M4 to M7), and replica R8 is assigned to modules 31a1 and 31a2 (M0 and M1). After four cycles, replica R4 is assigned to module 31a3 (M2), and replica R0 is assigned to module 31a4 (M3).

[0177] (S83a, S83b, S83c, S83d) DE calculations for the replicas assigned to each of the groups G0 to G3 are performed in parallel by the groups G0 to G3. When the number of modules M is 8, the DE calculations are performed by the ΔE calculation units 41c1 to 47cK shown in FIG.

[0178] (S84a, S84b, S84c, S84d) Flip determination for the replicas assigned to each of the groups G0 to G3 is performed in parallel by the groups G0 to G3. When the number of modules M is 8, the flip determination is performed by the selectors 40d to 47d shown in FIG. 4.

[0179] (S85a, S85b, S85c, S85d) Groups G0 to G3 determine in parallel whether a flip has occurred in the replica being processed by each group. Module control units 31b1 to 31bM of modules 31a1 to 31aM make the above determination based on flip bit information supplied from overall control unit 30.

[0180] If it is determined in the process of step S85a that a flip has occurred, the process of step S86a is performed. If it is determined in the process of step S85b that a flip has occurred, the process of step S86b is performed. If it is determined in the process of step S85c that a flip has occurred, the process of step S86c is performed. If it is determined in the process of step S85d that a flip has occurred, the process of step S86d is performed. If it is determined in the processes of steps S85a to S85d that a flip has not occurred, the process of step S88 is performed.

[0181] (S86a, S86b, S86c, S86d) In groups G0 to G3, each of modules 31a1 to 31aM reads all weighting coefficients related to flip bits from memory. If a flip bit occurs in all of groups G0 to G3, weighting coefficients related to the flip bits of all groups are read in each of modules 31a1 to 31aM (see FIG. 9).

[0182] (S87a, S87b, S87c, S87d) Each of the groups G0 to G3 uses the read weight coefficients to update h using the function shown in Fig. 5. If flip bits occur in all of the groups G0 to G3, the local fields corresponding to the bits for all groups are updated for each replica being processed in the groups G0 to G3.

[0183] (S88) Steps S82, S83a to S87a, S83b to S87b, S83c to S87c, and S83d to S87d are repeated until one trial has been performed for all bits in each replica. Once one trial has been performed for all bits in each replica, the loop is terminated and step S89 is performed.

[0184] (S89) The overall control unit 30 determines whether the search has ended. The overall control unit 30 determines that the search has ended when a predetermined search end condition is met. For example, the overall control unit 30 determines that the search has ended when the number of iterations reaches a predetermined number. If it is determined that the search has ended, the FPGA 28a ends the processing. If it is determined that the search has not ended, the processing from step S81 is repeated.

[0185] 10 to 15 are merely examples, and the order of the processes may be changed as appropriate. For example, the processes of steps S44 and S45 may be performed before the processes of steps S42 and S43.

[0186] According to the data processing device 20 of the second embodiment described above, the number of parallel trial bits P of the partial parallel trial is changed based on the search information indicating the search status. This makes it possible to set the number of parallel trial bits P according to the search status that reflects the characteristics of the problem, optimizes the amount of calculation required to change one bit, and improves the performance of solving large-scale problems.

[0187] In addition to optimizing the amount of calculation, by changing the number of parallel trial bits P as described above, when a bit is flipped, it is possible to adjust the period until that bit next becomes an update candidate. This makes it possible to avoid a situation where, when the state escapes from a local optimum due to the flipping of that bit, the bit is flipped again and the state is again constrained to a local optimum.

[0188] Furthermore, in the data processing device 20 of the second embodiment, n groups each including one or more modules perform parallel partial trials in parallel for n of the multiple replicas for each unit processing period (one step period). Among the groups performing the partial parallel trials, control is performed so that one group does not start processing the partial parallel trial for a replica until the update process (h update or state vector update) related to the number of parallel trial bits P for that replica is completed in the other groups. The data processing device 20 shifts the processing timing of the pipeline so that another group processes another replica until the update process for that replica is completed. This allows for effective use of computing resources while adhering to the principle of sequential processing of the MCMC method, even when the number of parallel trial bits P is variable, thereby improving the performance of solving relatively large-scale problems.

[0189] In the second embodiment, the number of groups is set to four as an example, but the number of groups may be a number other than four. The number of replicas may be a number other than 16. In addition, the number of bits handled by each module is set to K, but K may be a different value in each module.

[0190] Furthermore, the processing for each replica in the data processing device 20 may be executed by the FPGA 28a as in the above example, or may be executed by another computing unit such as the CPU 21 or the GPU 24. The computing units such as the FPGA 28a and the CPU 21 are examples of processing units in the data processing device 20. Furthermore, the storage unit that holds the multiple replicas may be realized by the memory 28b or a register as described above, or may be realized by the RAM 22. Furthermore, the accelerator card 28 can also be considered as an example of a "data processing device."

[0191] The information processing of the first embodiment may be realized by causing the processing unit 12 to execute a program. The information processing of the second embodiment may be realized by causing the CPU 21 to execute a program. The program can be recorded on a computer-readable recording medium 103.

[0192] For example, the program can be distributed by distributing the recording medium 103 on which it is recorded. Alternatively, the program may be stored in another computer and distributed via a network. A computer may store (install) a program recorded on the recording medium 103 or a program received from another computer in a storage device such as the RAM 22 or the HDD 23, and then read and execute the program from the storage device.

[0193] While one aspect of the program, data processing device, and data processing method of the present invention has been described above based on the embodiment, these are merely examples and the present invention is not limited to the above description. [Explanation of symbols]

[0194] 10 Data processing device 11 Storage section 12 Processing section

Claims

1. A program that causes a computer to execute a process of searching for a solution to a combinatorial optimization problem expressed as an energy function including a plurality of state variables, The computer, performing a search process for searching for the solution by determining in parallel whether or not to accept a change in the value of each of a plurality of first state variables selected from the plurality of state variables, and changing the value of any of the state variables for which it has been determined that the change in value is acceptable, while changing the selected plurality of first state variables; identifying the number of the plurality of first state variables to be selected based on search information indicating a search status of the search process or a search record of another combinatorial optimization problem, and repeating the search process; A program that executes a process.

2. The program according to claim 1 , wherein the search information includes a first cumulative value, during a first period, of the number of second state variables among the plurality of first state variables that are determined to accept the change in value.

3. calculating a first average value of the number of the second state variables during the first period based on the first cumulative value; If the first average value is greater than a first threshold, decreasing the number of the plurality of first state variables; increasing the number of the plurality of first state variables when the first average value is smaller than a second threshold value that is smaller than the first threshold value; The program according to claim 2 , which causes the computer to execute processing.

4. The program according to claim 1 , wherein the search information includes a second cumulative value of the number of third state variables whose values ​​have changed among the plurality of first state variables during a second period.

5. calculating an occurrence rate of the third state variable during the second time period based on the second cumulative value; If the occurrence rate is greater than a third threshold, decreasing the number of the plurality of first state variables; increasing the number of the plurality of first state variables when the occurrence rate is less than a fourth threshold that is less than the third threshold; The program according to claim 4, which causes the computer to execute processing.

6. The program according to claim 1 , wherein the search information includes a movement amount, expressed as a Hamming distance, of a state vector based on the plurality of state variables during a third period.

7. calculating a second average value of the movement amount during the third period based on the movement amount; If the second average value is greater than a fifth threshold and the minimum value of the energy function is not updated in the search process during the third period, decreasing the number of the plurality of first state variables; increasing the number of the plurality of first state variables when the second average value is smaller than a sixth threshold value that is smaller than the fifth threshold value and the minimum value is not updated; The program according to claim 6, which causes the computer to execute processing.

8. A data processing device that searches for a solution to a combinatorial optimization problem expressed by an energy function including a plurality of state variables, a processing unit that performs a search process for searching for the solution by determining in parallel whether or not to accept a change in the value of each of a plurality of first state variables selected from the plurality of state variables, changing the value of any of the state variables determined to accept the change in value while changing the plurality of selected first state variables, and identifies the number of the plurality of first state variables to be selected based on search information indicating a search status of the search process or a search record for another combinatorial optimization problem, and repeats the search process; a storage unit that stores the search information; A data processing device having:

9. the processing unit has M modules (M is an integer of 2 or more) that are grouped into n groups (n is an integer of 2 or more) each including one or more modules, and a selector; the n groups perform the determination on the plurality of first state variables in parallel for n replicas out of a plurality of replicas each indicating the plurality of state variables, for each unit processing period; the selector selects, in parallel, any of the state variables determined by the determination to accept the change in value for each of the n groups; the processing unit controls so that a group other than the first group among the n groups does not start processing of the first replica, which is one of the plurality of replicas, until an update process for changing the value of the state variable selected by the selector in the first group among the n groups is completed.

9. A data processing device according to claim 8.

10. A data processing method for searching for a solution to a combinatorial optimization problem expressed by an energy function including a plurality of state variables, the method comprising: performing a search process for searching for the solution by determining in parallel whether or not to accept a change in the value of each of a plurality of first state variables selected from the plurality of state variables, and changing the value of any of the state variables for which it has been determined that the change in value is acceptable, while changing the selected plurality of first state variables; identifying the number of the plurality of first state variables to be selected based on search information indicating a search status of the search process or a search record of another combinatorial optimization problem, and repeating the search process; Data processing methods.

Citation Information

Patent Citations

  • Optimization problem arithmetic program, optimization problem arithmetic method, and optimization problem arithmetic device

    JP2020046997A

  • Sampling device and sampling method

    JP2020204929A

  • Optimization device and control method of optimization device

    JP2021033341A

  • Optimization device, optimization method and control program of the optimization device

    JP2021131695A