WO2022249785A1 - 求解装置、求解方法およびプログラム - Google Patents
求解装置、求解方法およびプログラム Download PDFInfo
- Publication number
- WO2022249785A1 WO2022249785A1 PCT/JP2022/017778 JP2022017778W WO2022249785A1 WO 2022249785 A1 WO2022249785 A1 WO 2022249785A1 JP 2022017778 W JP2022017778 W JP 2022017778W WO 2022249785 A1 WO2022249785 A1 WO 2022249785A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- variable
- time
- elements
- unit
- updating
- Prior art date
Links
- 238000000034 method Methods 0.000 title claims description 57
- 238000005457 optimization Methods 0.000 claims abstract description 73
- 238000012937 correction Methods 0.000 claims abstract description 31
- 238000004364 calculation method Methods 0.000 claims description 79
- 230000008569 process Effects 0.000 claims description 44
- 230000010365 information processing Effects 0.000 claims description 22
- 230000008878 coupling Effects 0.000 claims description 6
- 238000010168 coupling process Methods 0.000 claims description 6
- 238000005859 coupling reaction Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000009471 action Effects 0.000 claims description 4
- 238000012545 processing Methods 0.000 abstract description 99
- 230000015654 memory Effects 0.000 description 32
- 238000007726 management method Methods 0.000 description 26
- 230000006870 function Effects 0.000 description 24
- 238000004891 communication Methods 0.000 description 18
- 238000010586 diagram Methods 0.000 description 14
- 230000005366 Ising model Effects 0.000 description 9
- 238000006243 chemical reaction Methods 0.000 description 4
- 238000009472 formulation Methods 0.000 description 4
- 239000000203 mixture Substances 0.000 description 4
- 230000001427 coherent effect Effects 0.000 description 3
- 230000003287 optical effect Effects 0.000 description 3
- 239000002245 particle Substances 0.000 description 3
- 238000002922 simulated annealing Methods 0.000 description 3
- 238000000137 annealing Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 239000000470 constituent Substances 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000005374 Kerr effect Effects 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000005401 electroluminescence Methods 0.000 description 1
- 230000005283 ground state Effects 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000005086 pumping Methods 0.000 description 1
- 238000012887 quadratic function Methods 0.000 description 1
- 238000013341 scale-up Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N99/00—Subject matter not provided for in other groups of this subclass
Definitions
- Embodiments of the present invention relate to a solution-seeking device, a solution-seeking method, and a program.
- Wm represents the upper limit.
- x A1 represents a decision variable that sets the investment ratio of stock A to 1%, for example.
- xA2 represents a decision variable that sets the investment ratio of stock A to 2%, for example.
- xSs represents a decision variable where the investment ratio of stock S is s%.
- ⁇ A1, ⁇ A2 and ⁇ Sm represent weights.
- the Tokyo Stock Exchange handles about 2000 stocks.
- the Tokyo Stock Exchange classifies these approximately 2,000 stocks into 33 types of industry and 3 types of scale. Therefore, the approximately 2,000 issues handled by the Tokyo Stock Exchange can be divided into a total of 99 groups of 33 industries x 3 scales.
- the number of decision variables included in the objective function is 200000.
- the number of Ising spins increases by about 50% in the case of the Tokyo Stock Exchange, for example, compared to when there are no constraints. Resulting in.
- the Ising problem is solved by a solver or the like, if the number of Ising spins increases, the amount of calculation and the calculation time will increase, and the cost for finding the solution will increase.
- the problem to be solved by the present invention is to provide a solution-finding device, a solution-finding method, and a program that can solve a 0-1 combinatorial optimization problem under one or more constraints represented by inequalities at a low computational cost. .
- a solution-seeking apparatus solves a 0-1 combinatorial optimization problem under one or more constraints represented by inequalities using a plurality of discrete variables included in the 0-1 combinatorial optimization problem.
- the solution-seeking device includes an updating unit and an output unit.
- the update unit alternately updates the first variable and the second variable sequentially from an initial time to an end time for each of the plurality of elements associated with the first variable and the second variable. do.
- the output unit outputs a solution of the 0-1 combinatorial optimization problem based on the first variable of each of the plurality of elements at the end time.
- the plurality of elements correspond to the plurality of discrete variables.
- Each of the first variable and the second variable is represented by a real number.
- the updating unit updates the first variable based on the second variable for each of the plurality of elements.
- the update unit updates the second variable based on the first variable for each of the plurality of elements.
- the updating unit performs the above-described Subtracting from the second variable of each of the plurality of elements, a correction value according to the component of the corresponding element in the distance from the boundary of the inequality to the position specified by the plurality of elements.
- FIG. 1 is a diagram illustrating the branching phenomenon of a simulated branching algorithm.
- FIG. 2 is a diagram showing the functional configuration of the solution-seeking device according to this embodiment.
- FIG. 3 is a flow chart showing a first example of the processing flow of the updating unit.
- FIG. 4 is a flowchart illustrating a second example of the processing flow of the updating unit.
- FIG. 5 is a flow chart showing the flow of processing of the first example of restriction processing.
- FIG. 6 is a flow chart showing the flow of processing of a second example of restriction processing.
- FIG. 7 is a flow chart showing the flow of a third example of restriction processing.
- FIG. 8 is a configuration diagram of an information processing system.
- FIG. 9 is a configuration diagram of the management server.
- FIG. 10 is a diagram illustrating data stored in a storage unit;
- FIG. 11 is a configuration diagram of a calculation server.
- FIG. 12 is a diagram showing data stored in a storage.
- An Ising machine is an example of a device used to solve the Ising problem.
- the Ising machine calculates the ground state energy of the Ising model.
- the Ising model has often been used mainly as a model for ferromagnets and phase transition phenomena.
- the Ising model has been increasingly used as a model for solving 0-1 combinatorial optimization problems. Equation (1) represents the energy of the Ising model.
- s i and s j represent spins.
- Spin is a binary variable that can take the value of either +1 or -1.
- s i represents the ith spin.
- s j represents the jth spin.
- i and j are integers of 1 or more and N or less.
- N represents the number of spins and is an integer of 2 or more.
- h i represents the local magnetic field acting on the i-th spin.
- J is a matrix of coupling coefficients representing the force acting between two spins.
- J is a real symmetric matrix with zeros on the diagonal.
- J ij represents the element of i row j column of J. That is, J ij is a coupling coefficient representing the force acting between the i-th spin and the j-th spin.
- the Ising machine uses the energy E Ising represented by Equation (1) as an objective function and calculates a solution that minimizes the energy E Ising .
- the Ising model solution (s 1 , s 2 , .
- the Ising model solution may not be the optimal solution but may be an approximate solution in which the energy E Ising is close to the minimum value. That is, the Ising problem may be a problem of calculating an approximate solution as well as an optimal solution.
- a 0-1 combinatorial optimization problem in which the objective function is a quadratic function of a discrete variable (bit) that takes a value of either 0 or 1 is called a 0-1 quadratic programming problem.
- Discrete variables (bits) are converted to s i by using the operation (1+s i )/2.
- the 0-1 quadratic programming problem is equivalent to the Ising problem represented by Equation (1). Therefore, the 0-1 quadratic programming problem can be converted into an Ising problem and the solution can be calculated by an Ising machine.
- Ising machines are hardware-implemented by, for example, quantum annealers, coherent Ising machines, and quantum bifurcation machines.
- Quantum annealers use superconducting circuits to achieve quantum annealing.
- a coherent Ising machine utilizes the oscillation phenomenon of a network formed of optical parametric oscillators.
- Quantum bifurcation machines exploit quantum mechanical bifurcation phenomena in networks of parametric oscillators with the Kerr effect.
- SA Simulated Annealing
- simulated annealing is a sequential update algorithm in which each variable is successively updated, it is difficult to speed up calculation processing by parallelization.
- Non-Patent Document 2 proposes a simulated branching algorithm as an algorithm for solving the 0-1 combinatorial optimization problem.
- the simulated bifurcation algorithm uses the Ising model and can solve a large-scale 0-1 combinatorial optimization problem at high speed by a digital computer.
- the simulated branching algorithm is a CPU (Central Processing Unit), a microprocessor, a GPU (Graphics Processing Unit), an FPGA (Field-Programmable Gate Array), an ASIC (Application Specific Integrated Circuit), or an electronic circuit such as a combination of these Circuits can also solve large-scale 0-1 combinatorial optimization problems at high speed.
- the simulated branching algorithm uses variables x i and y i each corresponding to N elements.
- the variable x i may be called the first variable, and the variable y i may be called the second variable.
- each of the N elements represents a virtual particle.
- the N elements correspond to the N spins of the Ising problem.
- the N elements thus correspond to the N discrete variables (bits) of the combinatorial optimization problem.
- Both the variables x i and the variables y i are continuous variables represented by real numbers.
- the variable x i represents the position of the i-th particle among the N particles.
- the variable yi represents the momentum of the ith particle.
- N is an integer of 2 or more.
- i represents an integer greater than or equal to 1 and less than or equal to N, and represents an index specifying each of the N elements.
- the simulated bifurcation algorithm numerically solves the simultaneous ordinary differential equations of equation (2) below for each of the N variables x i and y i .
- H is the Hamiltonian of formula (3) below.
- the coefficient D is a predetermined constant and corresponds to detuning.
- the coefficient p(t) corresponds to the pumping amplitude and increases monotonically with the number of updates during the calculation of the simulated branching algorithm.
- t is a variable representing time.
- the initial value of the coefficient p(t) may be set to zero.
- the coefficient K is a predetermined constant and corresponds to a positive Kerr coefficient. Note that the coefficient K may be 0.
- f i represents an external force and is represented by the following formula (4).
- Z i in equation (4) is an equation obtained by partially differentiating the expression in parentheses in equation (3) with variable x i .
- the expression in parentheses in equation (3) corresponds to the energy E Ising of the Ising model.
- c is a coefficient. c may be, for example, a constant that is predetermined before performing the calculation. Also, ⁇ (t) is a coefficient that increases with p(t).
- the simulated bifurcation algorithm then calculates the value of spin s i based on the sign of variable x i after increasing the value of p(t) from an initial value (eg, 0) to a predetermined value.
- the simulated bifurcation algorithm uses the symplectic Euler method to solve the differential equations given by equations (2), (3) and (4).
- t represents the time.
- ⁇ t represents unit time (time step, time interval).
- an electronic circuit such as a digital computer or an FPGA, based on the algorithm of equation (5) or equation (6), each of N variables x i and y i from the initial time unit time , and alternately updates the variables x i and y i . Then, an electronic circuit such as a digital computer or FPGA binarizes the values of N variables x i at the end time using a sign function, and outputs N spin values.
- equations (5) and (6) are expressed using time t and unit time ⁇ t in order to show correspondence with differential equations.
- the algorithm for computing equations (5) and (6) uses time t and unit time ⁇ t as explicit parameters. should not contain For example, if the unit time ⁇ t is 1, the algorithm for computing equations (5) and (6) may not include the unit time ⁇ t. For example, if we do not include time t as an explicit parameter, the algorithms that compute equations (5) and (6) operate with x i (t+ ⁇ t) as the updated value of x i (t). .
- FIG. 1 is a diagram showing a branching phenomenon of variables x i when an optimization problem is solved by a simulated branching algorithm.
- a bifurcation phenomenon occurs in which a system with only one stable motion state transitions to a system with two stable motion states as the parameters of the system change.
- the variables x i are concentrated in the vicinity of -1 or +1, but also spread into regions less than -1 or greater than +1.
- the improved algorithm numerically solves the system of ordinary differential equations of Eq. (2) using the Hamiltonian shown in Eq. (7) instead of the Hamiltonian of Eq. (3).
- G m represents a penalty term corresponding to the mth constraint of the one or more constraints.
- Each of the one or more constraints is represented by an inequality involving multiple discrete variables included in the objective function.
- the inequality is a first-order inequality that includes first-order terms in multiple discrete variables, or a second-order inequality that includes second-order terms in multiple discrete variables.
- the penalty term adds zero energy to the Hamiltonian if the corresponding inequality is satisfied, and adds positive energy to the Hamiltonian if the corresponding inequality is not satisfied.
- G m is a function that adds more energy to the Hamiltonian as the distance from the boundary of the constraint inequality to the position specified by the N elements increases. Note that the boundary of the inequality is the boundary between whether the inequality is satisfied or not. That is, the boundary of the inequality is represented by an equation in which the inequality sign is replaced by an equal sign.
- the improved algorithm solves the simultaneous ordinary differential equations by the symplectic Euler method using a discrete recurrence formula as shown in Equation (8) or Equation (9).
- the recurrence formulas shown in formulas (8) and (9) differ from the recurrence formulas shown in formulas (5) and (6) in the formula for calculating the second variable (y i (t+ ⁇ t)). Specifically, the recurrence formulas shown in formulas (8) and (9) add penalty terms ( G m ) is subtracted from the sum of the partial differential values of x i with respect to x i multiplied by ⁇ t.
- Each partial derivative of one or more inequalities corresponds to the component of the corresponding element at the distance from the boundary of the corresponding inequality to the position specified by the N elements.
- the penalty term becomes 0 when the corresponding inequality is satisfied. Therefore, the partial differential value of each of the one or more inequalities is 0 if N elements satisfy the inequalities.
- the simulated branching algorithm alternates between a plurality of first variables (x i ) and a plurality of second variables (y i ) while incrementing a variable t representing time every unit time ( ⁇ t). update to
- the simulated branching algorithm can determine, for each unit of time ( ⁇ t), whether an inequality is satisfied for each of the one or more constraints.
- the improved algorithm updates the plurality of second variables (y i ) every unit time ( ⁇ t), and then, for each of the one or more constraints, the first variable corresponding to each of the plurality of discrete variables It is determined whether or not the inequality into which (x i ) is substituted is satisfied. Then, if the improved algorithm does not satisfy the inequality obtained by substituting the first variable corresponding to each of the plurality of discrete variables for each of the one or more constraints for each unit time ( ⁇ t), then the plurality of second variables From each of (y i ), subtract the corresponding correction value.
- the plurality of second variables From each of (y i ), subtract a correction value according to the component of the corresponding element in the distance from the boundary of the inequality to the position specified by the multiple elements.
- the improved algorithm can generate a plurality of first variables (x i ) and a plurality of second variables (y i ) for each unit time ( ⁇ t) by a plurality of A first variable (x i ) and a plurality of second variables (y i ) can be corrected.
- Such an improved algorithm solves a 0-1 combinatorial optimization problem under one or more constraints expressed by inequalities, with the same number of variables as when no constraints are set. , the solution can be obtained. Therefore, such an improved algorithm can easily solve a 0-1 combinatorial optimization problem under one or more constraints expressed by inequalities without providing a slug variable.
- One of the one or more constraints is expressed by a linear inequality.
- the improved algorithm solves a 0-1 combinatorial optimization problem under the constraints of M linear inequalities, where M is an integer greater than or equal to 1.
- the m-th constraint among the one or more constraints is represented by Equation (10).
- m is an integer of 1 or more and M or less.
- x i is the i-th discrete variable.
- x j is the jth discrete variable.
- ⁇ m,i is the coefficient multiplied by the i th discrete variable.
- Wm is a predetermined constant.
- a m is a predetermined coefficient for the mth constraint.
- k is an integer of 1 or more.
- k may be two.
- S m is a value obtained by substituting the first variable (x i ) corresponding to each of the plurality of discrete variables in the function except for the constant term (W m ) in the first-order inequality. Specifically, Sm is expressed as in Equation (12) when the recurrence formula of Equation (8) is applied. Moreover, Sm is represented by Formula (13) when applying the recurrence formula of Formula (9).
- Equation (14) the partial differential value is expressed as shown in Equation (14).
- the improved algorithm updates the plurality of second variables (y i ) every unit time ( ⁇ t), and then updates the first It is determined whether or not the linear inequality substituting the variable (x i ) is satisfied. Then, if the corresponding linear inequality is not satisfied for each unit time ( ⁇ t), the improved algorithm converts each of the plurality of second variables (y i ) into the partial differential value represented by Equation (14) in the unit time Subtract the correction value multiplied by ⁇ t representing .
- An improved algorithm that performs such a process can easily solve the 0-1 combinatorial optimization problem under the constraint expressed by the linear inequality with the same number of variables as when no constraint is set. .
- Equation (15) when the constraint is a linear inequality. That is, the improved algorithm converts each of the plurality of second variables (y i ) into the partial differential value represented by Equation (15) for each unit time ( ⁇ t) if the corresponding linear inequality is not satisfied. A correction value multiplied by ⁇ t representing , may be subtracted.
- B m is a value greater than or equal to 0 and less than or equal to 1, and is a predetermined coefficient for the m-th constraint.
- E m is a value obtained by substituting the corresponding second variable (y i ) for each of a plurality of discrete variables in the function excluding the constant term (W m ) in the first-order inequality. Specifically, E m is expressed as in Equation (16) when the recurrence formula of Equation (8) or Equation (9) is applied.
- E m represents the slope of the direction specified by the elements with respect to the boundary of the linear inequality. Therefore, when the partial differential value of equation (15) is used, the correction value is the component of the corresponding element in the distance from the boundary of the inequality to the position specified by the multiple elements and the direction specified by the multiple elements. , and the component of the element corresponding to the gradient with respect to the boundary of the first-order inequality. This allows the improved algorithm to stabilize the correction value even when the component with distance becomes unstable due to large Am .
- Equation (17) (Penalty term and partial differential value of constraint is quadratic inequality)
- Equation (17) One of the one or more constraints may be expressed by a quadratic inequality.
- the constraint is represented by Equation (17).
- x i is the i-th discrete variable.
- x j is the jth discrete variable.
- Q i,j is the value of i-th row, j-th column contained in the N-th row, N-column positive semidefinite matrix.
- q is a predetermined constant.
- A is a predetermined coefficient for the constraint condition by the quadratic inequality.
- k is an integer of 1 or more.
- k may be two.
- S is a value obtained by substituting the corresponding first variables (x i , x j ) for each of the plurality of discrete variables in the function excluding the constant term (q) in the quadratic inequality.
- S is represented by formula (19).
- S is represented like Formula (20), when applying the recurrence formula of Formula (9).
- the partial differential value is expressed as shown in Equation (21) when applying the recurrence formula of Equation (8). Also, the partial differential value is expressed as shown in Equation (22) when the recurrence formula of Equation (9) is applied.
- the improved algorithm updates the plurality of second variables (y i ) every unit time ( ⁇ t), and then updates the plurality of second variables (y i ), and then the second variables corresponding to each of the plurality of discrete variables. It is determined whether or not the quadratic inequality substituting one variable (x i , x j ) is satisfied. Then, the improved algorithm is represented by equation (21) or equation (22) from each of the plurality of second variables (y i ) if the corresponding quadratic inequality is not satisfied for each unit time ( ⁇ t) A correction value obtained by multiplying the partial differential value by ⁇ t representing unit time is subtracted.
- the improved algorithm can easily solve the 0-1 combinatorial optimization problem under the constraint expressed by the quadratic inequality.
- the improved algorithm can solve a 0-1 combinatorial optimization problem even if one or more of the constraints include both linear and quadratic inequality constraints. can.
- the improved algorithm may also include only linear inequality constraints in one or more of the constraints.
- the improved algorithm may also solve the 0-1 combinatorial optimization problem under only one quadratic inequality constraint.
- FIG. 2 is a diagram showing the functional configuration of the solution-seeking device 10 according to this embodiment.
- the solving device 10 uses an improved algorithm to solve the 0-1 combinatorial optimization problem under one or more constraints represented by inequalities.
- the solution-finding device 10 is an information processing device such as a computer, a computer system configured by a plurality of computers or servers communicating with each other via a network, or a PC in which a plurality of computers cooperate to execute information processing. It is implemented by a cluster or the like. Also, the solution-seeking device 10 is realized by an electronic circuit such as a CPU, a microprocessor, a GPU, an FPGA or an ASIC, or a circuit of a combination thereof.
- the solution-seeking device 10 includes an input unit 12, an update unit 14, and an output unit 16 as functional configurations.
- the input unit 12 receives information (eg, N, J, h) for defining the objective function of the 0-1 combinatorial optimization problem, and information (eg, D , c, ⁇ t, T, p(t), ⁇ (t)) from an external device. Furthermore, the input unit 12 receives information (for example, ⁇ m,i , W m , Q i,j , q) for defining one or more constraint conditions, and coefficients necessary for executing constraint processing. Information representing (eg, A m , k, B m , A) is received from an external device. The input unit 12 then provides the received information to the updating unit 14 .
- information eg, N, J, h
- information eg, D , c, ⁇ t, T, p(t), ⁇ (t)
- information for example, ⁇ m,i , W m , Q i,j , q
- Information representing eg, A m ,
- each of the first variable (x i ) and the second variable (y i ) is represented by a real number.
- the updating unit 14 updates the first variable (x i ) based on the second variable (y i ) for each of the plurality of elements.
- the update unit 14 updates the second variable (y i ) based on the first variable (x i ) for each of the plurality of elements.
- the updating unit 14 updates the second variable (y i ) after updating the first variable (x i ) for each of the plurality of elements.
- the updating unit 14 may update the first variable (x i ) after updating the second variable (y i ) for each of the plurality of elements.
- the updating unit 14 performs constraint processing so that the solution of the 0-1 combinatorial optimization problem satisfies one or more constraint conditions after updating the second variable (y i ). to run. For each of the one or more constraints, the updating unit 14 updates the second variables ( y i ), subtract the correction value (r m,i ) according to the component of the corresponding element at the distance from the boundary of the corresponding inequality to the position specified by the plurality of elements.
- the updating unit 14 updates the second variable Subtract the correction value (r m,i ) from the second variable (y i ) after updating (y i ).
- the updating unit 14 updates the second variable Subtract the correction value (r m, i ) from the second variable (y i ) after updating (y i ) and before updating the first variable (x i ).
- FIG. 3 is a flow chart showing a first example of the processing flow of the updating unit 14. As shown in FIG. The updating unit 14 executes processing according to the flow shown in FIG. 3, for example.
- the updating unit 14 sets information about the constraint. Specifically, the updating unit 14 sets information for defining one or more constraint conditions and information representing coefficients necessary for executing constraint processing. For example, when the constraint condition is represented by a linear inequality, the updating unit 14 updates ⁇ m,i , which is a coefficient to be multiplied by the variable, W m , which is a constant term, and Set A m (or A m and B m ) and k.
- the update unit 14 updates Qi ,j , which is a coefficient to be multiplied by the variable, q, which is a constant term, and a coefficient necessary for executing the constraint process.
- Qi ,j which is a coefficient to be multiplied by the variable, q, which is a constant term, and a coefficient necessary for executing the constraint process.
- Set A and k such that Note that the update unit 14 may set A m , B m , k, and A, which are coefficients necessary for executing the constraint processing, according to the parameters received from the input unit 12, or may be determined in advance. You can also set a value that is fixed and cannot be changed.
- the updating unit 14 initializes variables. Specifically, the updating unit 14 initializes a variable t representing time to an initial time (for example, 0). Furthermore, the updating unit 14 updates each of the N first variables (x 1 (t) to x N (t)) and each of the N second variables (y 1 (t) to y N (t)) is substituted with an initial value received from the user, a predetermined fixed value, or a random number.
- a variable t representing time to an initial time (for example, 0).
- the updating unit 14 updates each of the N first variables (x 1 (t) to x N (t)) and each of the N second variables (y 1 (t) to y N (t)) is substituted with an initial value received from the user, a predetermined fixed value, or a random number.
- the update unit 14 repeats the loop processing between S104 and S116 until t becomes greater than T.
- the update unit 14 converts the N first variables (x 1 (t+ ⁇ t) to x N (t+ ⁇ t)) at the target time (t+ ⁇ t) to the N first variables at the immediately preceding time (t). It is calculated based on variables (x 1 (t) to x N (t)) and N second variables (y 1 (t) to y N (t)) at the immediately preceding time (t).
- the updating unit 14 converts N second variables (y 1 (t+ ⁇ t) to y N (t+ ⁇ t)) at the target time (t+ ⁇ t) to N second variables at the target time (t+ ⁇ t). It is calculated based on the first variables (x 1 (t+ ⁇ t) to x N (t+ ⁇ t)) and N second variables (y 1 (t) to y N (t)) at the immediately preceding time (t).
- i is an integer from 1 to N and is an index representing a processing target among N elements.
- Each of the N elements is associated with a first variable (x i (t)) and a second variable (y i (t)).
- the update unit 14 performs processing with the i-th element of the N elements as the target element.
- the updating unit 14 converts the first variable (x i (t+ ⁇ t)) at the target time (t+ ⁇ t) of the target element to the first variable (x i (t)) at the immediately preceding time (t) of the target element, It is calculated by adding the value obtained by multiplying the second variable (y i (t)) at the time (t) immediately before the target element, a predetermined constant (D), and the unit time ( ⁇ t). Specifically, the updating unit 14 calculates Equation (23).
- the updating unit 14 converts the first variable (x i (t+ ⁇ t)) at the target time (t+ ⁇ t) of the target element to the first variable (x i (t)) and the second variable (y i (t)) at the previous time (t) of the target element.
- the update unit 14 updates the first variables (x 1 (t+ ⁇ t) to x N (t+ ⁇ t)) at the target time (t+ ⁇ t) of each of the N elements, and the relationship between the target element and each of the N elements.
- An update value (z i (t+ ⁇ t)) is calculated based on an action coefficient predetermined by a 0-1 combinatorial optimization problem for each set.
- the action coefficients are the coupling coefficients contained in J and the local magnetic field coefficients contained in h.
- the updating unit 14 calculates Equation (24).
- the update unit 14 calculates the external force (f i (t+ ⁇ t)) by multiplying the update value (z i (t+ ⁇ t)) by the coefficient (c) and ⁇ 1. Specifically, the updating unit 14 calculates Equation (25).
- the update unit 14 multiplies the value determined based on p(t+ ⁇ t), which is a function that increases with the passage of time, by the first variable (x 1 (t+ ⁇ t)) at the target time (t+ ⁇ t) of the target element. Then, the time evolution value (g i (t+ ⁇ t)) is calculated. Specifically, the updating unit 14 calculates Equation (26).
- the updating unit 14 converts the second variable (y i (t+ ⁇ )) at the target time (t+ ⁇ t) of the target element to the second variable (y i (t) ), a value obtained by multiplying the sum of the time evolution value (g i (t+ ⁇ t)) and the external force (f i (t+ ⁇ t)) by the unit time ( ⁇ t). Specifically, the updating unit 14 calculates Equation (27).
- the updating unit 14 executes the above-described loop processing between S108 and S113 N times, thereby updating the second variable (y i (t+ ⁇ t)) at the target time (t+ ⁇ t) for each of the N elements.
- N first variables x 1 (t + ⁇ t) to x N (t + ⁇ t)) at the target time (t + ⁇ t)
- the second variable y i (t)) at the time (t) immediately before the target element Update based on.
- the update unit 14 executes constraint processing based on one or more constraint conditions. Details of the restriction process will be described later with reference to the flowcharts of FIGS. 5 to 7. FIG. After completing the process of S114, the updating unit 14 advances the process to S115.
- the update unit 14 adds the unit time ( ⁇ t) to each of the previous time (t) and the target time (t+ ⁇ t) to update the previous time (t) and the target time (t+ ⁇ t).
- the update unit 14 repeats the processing from S105 to S115 until t exceeds the end time (T). Then, when t becomes greater than the end time (T), the updating unit 14 ends this flow.
- the update unit 14 executes the calculation according to the improved algorithm, and the N first variables (x 1 (t) to x N (t)) and N second variables (y 1 (t) to y N (t)) can be calculated.
- FIG. 4 is a flowchart showing a second example of the processing flow of the updating unit 14.
- the updating unit 14 may execute the process according to the flow shown in FIG. 4 instead of the flow shown in FIG.
- the update unit 14 performs the same processing as S101, S102 and S103 of the first example shown in FIG.
- the update unit 14 repeats the loop processing between S204 and S216 until t becomes greater than T.
- the update unit 14 converts N second variables (y 1 (t+ ⁇ t) to y N (t+ ⁇ t)) at the target time (t+ ⁇ t) to N first variables at the previous time (t). It is calculated based on variables (x 1 (t) to x N (t)) and N second variables (y 1 (t) to y N (t)) at the immediately preceding time (t).
- the update unit 14 updates the N first variables (x 1 (t+ ⁇ t) to x N (t+ ⁇ t)) at the target time (t+ ⁇ t) to the N first variables at the immediately preceding time (t). It is calculated based on the first variables (x 1 (t) to x N (t)) and N second variables (y 1 (t+ ⁇ t) to y N (t+ ⁇ t)) at the target time (t+ ⁇ t).
- the updating unit 14 performs processing with the i-th element among the N elements as the target element.
- the updating unit 14 updates the first variables (x 1 (t) to x N (t)) at the time (t) immediately before each of the N elements, and the target element and each of the N elements.
- An update value (z i (t)) is calculated based on an action coefficient predetermined by a 0-1 combinatorial optimization problem for each set. Specifically, the updating unit 14 calculates Equation (28).
- the update unit 14 calculates the external force (f i (t)) by multiplying the updated value (z i (t)) by the coefficient (c) and ⁇ 1. Specifically, the updating unit 14 calculates Equation (29).
- the updating unit 14 multiplies the value determined based on p(t), which is a function that increases over time, by the first variable (x 1 (t)) at the time (t) immediately before the target element. Then, the time evolution value (g i (t)) is calculated. Specifically, the updating unit 14 calculates Equation (30).
- the updating unit 14 converts the second variable (y i (t+ ⁇ )) at the target time (t+ ⁇ t) of the target element to the second variable (y i (t) ), a value obtained by multiplying the sum of the time evolution value (g i (t)) and the external force (f i (t)) by the unit time ( ⁇ t). Specifically, the updating unit 14 calculates Equation (31).
- the update unit 14 executes the above-described loop processing between S205 and S210 N times, thereby updating the second variable (y i (t+ ⁇ t)) at the target time (t+ ⁇ t) for each of the N elements.
- the update unit 14 executes constraint processing based on one or more constraint conditions. Details of the restriction process will be described later with reference to the flowcharts of FIGS. 5 to 7. FIG. After completing the process of S211, the update unit 14 advances the process to S212.
- the updating unit 14 converts the first variable (x i (t+ ⁇ t)) at the target time (t+ ⁇ t) of the target element to the first variable (x i (t)) at the immediately previous time (t) of the target element, It is calculated by adding a value obtained by multiplying a second variable (y i (t+ ⁇ t)), a predetermined constant (D), and a unit time ( ⁇ t) at the object time (t+ ⁇ t) of the object element. Specifically, the updating unit 14 calculates Equation (32).
- the updating unit 14 converts the first variable (x i (t+ ⁇ t)) at the target time (t+ ⁇ t) of the target element to the first variable (x i (t)) and the second variable (y i (t)) at the previous time (t) of the target element.
- the update unit 14 adds the unit time ( ⁇ t) to each of the previous time (t) and the target time (t+ ⁇ t) to update the previous time (t) and the target time (t+ ⁇ t).
- the updating unit 14 repeats the processing from S212 to S215 until t exceeds the end time (T). Then, when t becomes greater than the end time (T), the updating unit 14 ends this flow.
- the updating unit 14 executes calculations according to the improved algorithm to obtain N first variables (x 1 (t) to x N (t)) and N second variables (y 1 (t) to y N (t)) can be calculated.
- FIG. 5 is a flow chart showing the flow of processing of the first example of restriction processing.
- the updating unit 14 performs the processing in S114 of FIG. 3 and S211 of FIG. 4 according to the flow shown in FIG.
- m is an index representing a primary inequality to be processed among the M primary inequalities.
- the updating unit 14 calculates the constraint value (S m ). Specifically, in the constraint processing of S114 of FIG. 3, the updating unit 14 performs the calculation of Expression (33). In addition, in the constraint processing of S211 in FIG. 4, the updating unit 14 performs the calculation of Expression (34).
- the control value S m is a value obtained by substituting the corresponding first variable (x i ) for each of the plurality of discrete variables included in the function excluding the constant term (W m ) in the corresponding linear inequality. is.
- the update unit 14 determines whether or not the constraint value (S m ) is equal to or less than the constant term (W m ) in the corresponding linear inequality. That is, the updating unit 14 determines whether or not the plurality of first variables (x i ) satisfy corresponding linear inequalities.
- the updating unit 14 performs processing with the i-th element among the N elements as the target element.
- the updating unit 14 calculates a correction value (r m, i ) is calculated. Specifically, the update unit 14 performs the calculation of Expression (35).
- the updating unit 14 subtracts the correction value (r m,i ) from the second variable (y i ) of the target element. Specifically, the updating unit 14 performs the calculation of Equation (36).
- the update unit 14 can calculate the second A correction value (r m,i ) depending on the component of the target element at the distance from the boundary of the linear inequality to the position specified by the multiple elements can be subtracted from the variable (y i ).
- the update unit 14 changes in the direction in which the plurality of first variables (x i ) satisfy the corresponding linear inequalities.
- a plurality of second variables (y i ) can be corrected.
- FIG. 6 is a flow chart showing the flow of processing of a second example of restriction processing.
- the updating unit 14 performs the processing in S114 of FIG. 3 and S211 of FIG. 4 according to the flow shown in FIG. good too.
- the updating unit 14 calculates the constraint value (S m ).
- S m the constraint value
- the updating unit 14 determines whether or not the constraint value (S m ) is equal to or less than the constant term (W m ) in the corresponding linear inequality.
- the processing of S403 is the same as the processing of S303 in FIG. If the constraint value (S m ) is equal to or less than the constant term (W m ) in the corresponding linear inequality (Yes in S403), the updating unit 14 advances the process to S409. If the constraint value (S m ) is not equal to or less than the constant term (W m ) in the corresponding linear inequality (No in S403), the updating unit 14 advances the process to S404.
- the update unit 14 calculates the slope (E m ) with respect to the boundary of the corresponding linear inequality in the direction specified by the multiple elements. Specifically, the updating unit 14 performs the calculation of Expression (37).
- the updating unit 14 performs processing with the i-th element of the N elements as the target element.
- the update unit 14 determines the component of the target element at the distance from the boundary of the linear inequality to the position specified by the multiple elements and the A correction value (r m,i ) is calculated according to the component of the target element of the inclination with respect to the boundary of the direction linear inequality. Specifically, the updating unit 14 performs the calculation of Expression (38).
- the update unit 14 subtracts the correction value (r m,i ) from the second variable (y i ) of the target element.
- the processing of S407 is the same as the processing of S306 in FIG.
- the updating unit 14 can calculate the second From the variable (y i ), the component of the target element at the distance from the boundary of the linear inequality to the position specified by the plurality of elements, and the slope of the target element with respect to the boundary of the primary inequality in the direction specified by the plurality of elements A component-dependent correction value (r m,i ) can be subtracted.
- the update unit 14 adds to the correction value (r m,i ) a value corresponding to the component of the target element of the gradient with respect to the boundary of the linear inequality in the direction specified by the plurality of elements, so that the unit time ( ⁇ t ), it is possible to suppress unstable fluctuations of the plurality of first variables (x i ) and the plurality of second variables (y i ).
- the updating unit 14 stably changes in the direction in which the plurality of first variables (x i ) satisfy the corresponding linear inequalities.
- a plurality of second variables (y i ) corresponding to a plurality of elements can be corrected such that:
- FIG. 7 is a flow chart showing the flow of a third example of restriction processing. If a quadratic inequality is set for one of the one or more constraints, the update unit 14 performs the quadratic inequality instead of the processing of S302 to S307 in FIG. 5 or S402 to S408 in FIG. Then, the processing of S501 to S506 shown in FIG. 7 is executed.
- the update unit 14 calculates a constraint value (S). Specifically, in the constraint processing of S114 of FIG. 3, the updating unit 14 performs the calculation of Expression (39). In addition, in the constraint processing of S211 of FIG. 4, the updating unit 14 executes the calculation of Expression (40).
- the control value S is a value obtained by substituting the corresponding first variable (x i ) for each of the plurality of discrete variables included in the function excluding the constant term (q) in the quadratic inequality. .
- the update unit 14 determines whether or not the constraint value (S) is equal to or less than the constant term (q) in the quadratic inequality. That is, the updating unit 14 determines whether or not the plurality of first variables (x i ) satisfy the quadratic inequality.
- the updating unit 14 performs processing with the i-th element among the N elements as the target element.
- the updating unit 14 calculates the correction value (r i ) according to the component of the target element at the distance from the boundary of the quadratic inequality to the position specified by the plurality of elements. Calculate Specifically, in the constraint processing of S114 of FIG. 3, the updating unit 14 performs the calculation of Equation (41). In addition, in the constraint processing of S211 of FIG. 4, the updating unit 14 performs the calculation of Expression (42).
- the update unit 14 subtracts the correction value (r i ) from the second variable (y i ) of the target element. Specifically, the updating unit 14 performs the calculation of Expression (43).
- the updating unit 14 can obtain the second variables (y From i ), a correction value (r i ) depending on the component of the object element of the distance from the boundary of the quadratic inequality to the position specified by the multiple elements can be subtracted.
- the updating unit 14 changes in the direction in which the plurality of first variables (x i ) satisfy the corresponding quadratic inequalities.
- multiple second variables (y i ) can be corrected.
- the solution-seeking device 10 solves the 0-1 combinatorial optimization problem using the simulated branching algorithm.
- the solution-finding device 10 according to the present embodiment does not satisfy the inequality in which the first variable corresponding to each of the plurality of discrete variables is substituted for each of one or more constraints. , subtract the correction value (r i ) according to the component of the corresponding element in the distance from the boundary of the inequality to the position specified by the plurality of elements from the second variable (y i ) of each of the plurality of elements .
- the solution-seeking apparatus 10 can solve the 0-1 combinatorial optimization problem under one or more constraints represented by inequalities, even when no constraints are set.
- a solution can be obtained with the same number of variables. Therefore, according to the solving device 10 according to the present embodiment, it is possible to easily solve a 0-1 combinatorial optimization problem under one or more constraints represented by inequalities without providing slug variables. .
- the solution-seeking device 10 according to the present embodiment can reduce the amount of memory for storing intermediate progress of variables, etc., and can reduce the amount of calculation and the calculation time.
- the solution-finding device 10 when executing the constraint process of the linear inequality, the solution-finding device 10 generates M linear inequalities, and ⁇ m,i , which is a coefficient to be multiplied by the i-th discrete variable, and the constant of the linear inequality You may call some W m .
- the solution-seeking apparatus 10 can simultaneously read ⁇ m, i and W m by specifying the column number and row number once, and can easily execute processing.
- FIG. 8 is a diagram showing the configuration of the information processing system 100. As shown in FIG. The refinement algorithm can be executed by an information processing system 100 as shown in FIG. 8, for example. As a result, the information processing system 100 can solve large-scale combinatorial optimization problems at high speed through parallel processing.
- the information processing system 100 includes a management server 101, a network 102, a plurality of calculation servers 103 (103a to 103c) (information processing apparatuses), a plurality of cables 104 (104a to 104c), and a switch 105. Also, the terminal device 106 can communicate with the information processing system 100 .
- the management server 101, the plurality of calculation servers 103 (103a to 103c), and the terminal device 106 perform data communication with each other via the network 102.
- FIG. Network 102 is, for example, the Internet in which multiple computer networks are interconnected. Network 102 may be wired, wireless, or a combination thereof as communication media.
- a plurality of calculation servers 103 are connected to the switch 105 via cables 104 (104a to 104c) respectively.
- a plurality of cables 104 (104a-104c) and switches 105 form the interconnects between the computation servers 103 and the computation servers 103 .
- Each of the plurality of calculation servers 103 (103a to 103c) can also perform data communication with each other via interconnects.
- the switch 105 is, for example, an Infiniband switch. Cables 104a-104c are, for example, Infiniband cables. However, switch 105 and cable 104 may be a wired LAN switch/cable.
- the communication standard and communication protocol used by the cable 104 and the switch 105 are not particularly limited.
- the terminal device 106 is, for example, a notebook PC, a desktop PC, a smart phone, a tablet, or an in-vehicle terminal device.
- each of the plurality of calculation servers 103 (103a to 103c) and/or the processors of the calculation servers 103 (103a to 103c) perform part of the calculation processing in solving the combinatorial optimization problem using the improved algorithm. , or the same calculation process for different variables may be executed in parallel.
- the management server 101 converts the combinatorial optimization problem input by the user into a format that can be processed by each calculation server 103, and causes each calculation server 103 to execute it. Then, the management server 101 acquires the calculation results from each calculation server 103 and converts the aggregated calculation results into a solution of the combinatorial optimization problem.
- FIG. 8 shows three calculation servers 103 (103a to 103c).
- the number of calculation servers 103 included in the information processing system 100 is not limited to three.
- the information processing system 100 may be one, two, or four or more.
- the information processing system 100 may use some of the plurality of calculation servers 103 included in the process of solving a combinatorial optimization problem.
- Calculation server 103 may be any type of information processing device.
- the calculation server 103 may be a server installed in a data center or a desktop PC installed in an office.
- the calculation server 103 may be a plurality of types of computers installed at different locations.
- computation server 103 may be a general-purpose computer, dedicated electronic circuitry, or a combination thereof.
- FIG. 9 is a diagram showing the configuration of the management server 101.
- the management server 101 is, for example, a computer including a central processing unit (CPU) and memory.
- the management server 101 includes a processor 110 , a storage unit 114 , a communication circuit 115 , an input circuit 116 and an output circuit 117 .
- Processor 110 , storage unit 114 , communication circuit 115 , input circuit 116 and output circuit 117 are connected to each other via bus 120 .
- the processor 110 includes a management unit 111, a conversion unit 112, and a control unit 113 as an internal functional configuration.
- the processor 110 is an electronic circuit that performs calculations and controls the management server 101 .
- Processor 110 may be, for example, a CPU, microprocessor, ASIC, FPGA, PLD, or a combination thereof.
- the management unit 111 provides an interface for operating the management server 101 via the terminal device 106 of the user.
- the interface provided by the management unit 111 is, for example, API, CLI, web page, or the like.
- the user inputs information on the combinatorial optimization problem via the management unit 111, and browses and/or downloads the calculated solution of the combinatorial optimization problem.
- the conversion unit 112 inputs parameters relating to the combinatorial optimization problem and converts the input parameters into a format that can be processed by each calculation server 103 .
- the control unit 113 transmits a control command to each calculation server 103 .
- the conversion unit 112 aggregates a plurality of calculation results, converts them into a solution of the combinatorial optimization problem, and outputs the solution of the combinatorial optimization problem.
- the storage unit 114 stores various data including programs of the management server 101, data necessary for executing the programs, and data generated by the programs.
- a program includes both an OS and an application.
- Storage unit 114 may be volatile memory, non-volatile memory, or a combination thereof. Volatile memory is, for example, DRAM or SRAM.
- the nonvolatile memory is NAND flash memory, NOR flash memory, ReRAM, MRAM, or the like.
- the storage unit 114 may be a hard disk, an optical disk, a magnetic tape, an external storage device, or the like.
- the communication circuit 115 transmits and receives data to and from each device connected to the network 102 .
- the communication circuit 115 is, for example, a wired LAN NIC (Network Interface Card). However, the communication circuit 115 may be another type of communication circuit such as a wireless LAN.
- the input circuit 116 implements data input to the management server 101 .
- the input circuit 116 includes, for example, USB, PCI-Express, etc. as external ports.
- the operating device 118 is connected to the input circuit 116 .
- the operation device 118 is a device for the user to input information to the management server 101 .
- the operation device 118 is, for example, a keyboard, mouse, touch panel, voice recognition device, etc., but is not limited to these.
- the output circuit 117 implements data output from the management server 101 .
- the output circuit 117 includes HDMI (registered trademark), DisplayPort, or the like as an external port.
- display device 119 is connected to output circuit 117 .
- the display device 119 is, for example, an LCD (liquid crystal display), an organic EL (organic electroluminescence) display, a projector, or the like, but is not limited to these.
- the administrator can use the operation device 118 and the display device 119 to perform maintenance on the management server 101.
- the operation device 118 and the display device 119 may be incorporated in the management server 101 .
- the operation device 118 and the display device 119 do not have to be connected to the management server 101 .
- the administrator may maintain the management server 101 using a terminal device capable of communicating with the network 102 .
- FIG. 10 is a diagram showing data stored in the storage unit 114 of the management server 101.
- the storage unit 114 stores, for example, question data 114A, calculation data 114B, a management program 114C, a conversion program 114D, and a control program 114E.
- problem data 114A includes data for combinatorial optimization problems.
- the calculation data 114B includes calculation results collected from each calculation server 103 .
- the management program 114C is a program that implements the functions of the management unit 111 .
- the conversion program 114D is a program that implements the functions of the conversion unit 112 .
- the control program 114E is a program that implements the functions of the control unit 113 .
- FIG. 11 is a diagram showing the configuration of the calculation server 103a.
- the other calculation server 103 may have the same configuration as the calculation server 103a, or may have a different configuration from the calculation server 103a.
- the calculation server 103a includes, for example, a communication circuit 131, a shared memory 132, processors 133a to 133d, a storage 134, and a host bus adapter 135.
- Communication circuit 131 , shared memory 132 , processors 133 a to 133 d , storage 134 and host bus adapter 135 are connected to each other via bus 136 .
- the communication circuit 131 transmits and receives data to and from each device connected to the network 102 .
- the communication circuit 131 is, for example, a wired LAN NIC (Network Interface Card). However, the communication circuit 131 may be another type of communication circuit such as a wireless LAN.
- the shared memory 132 is memory accessible by the processors 133a-133d.
- the shared memory 132 is, for example, volatile memory such as DRAM and SRAM.
- Shared memory 132 may include other types of memory, such as non-volatile memory.
- Processors 133 a - 133 d share data via shared memory 132 . It should be noted that the shared memory 132 may not be composed of all the memory of the calculation server 103a. For example, some memory of computation server 103a may be local memory that is accessible only by any of processors 133a-133d.
- the processors 133a-133d are electronic circuits that perform calculation processing.
- the processors 133a-133d may be, for example, CPUs, GPUs, FPGAs, ASICs, or combinations thereof.
- the processors 133a-133d may be CPU cores or CPU threads. When the processors 133a to 133d are CPUs, the number of sockets provided in the calculation server 103a does not matter.
- the processors 133a-133d may also be connected to other components of the computation server 103a via buses such as PCI express.
- the calculation server 103a comprises four processors 133a-133d.
- the number of processors included in one calculation server 103a is not limited to four.
- the storage 134 stores various data including programs of the calculation server 103a, data necessary for executing the programs, and data generated by the programs.
- the program includes both OS and applications.
- Storage 134 may be volatile memory and non-volatile memory, or a combination thereof. Volatile memory is, for example, DRAM or SRAM. Non-volatile memory is, for example, NAND flash memory, NOR flash memory, ReRAM, MRAM, or the like. Storage 134 may also include hard disks, optical disks, magnetic tapes, or external storage devices.
- the host bus adapter 135 implements data communication with other calculation servers 103 .
- Host bus adapter 135 is connected to switch 105 via cable 104a.
- the host bus adapter 135 is, for example, an HCA (Host Channel Adapter).
- HCA Hyper Channel Adapter
- FIG. 12 is a diagram showing data stored in the storage 134.
- FIG. Storage 134 stores, for example, calculation data 134A, calculation program 134B, and control program 134C.
- the calculation data 134A includes data during calculation or calculation results of the calculation server 103a. At least a portion of computational data 134A may be stored in different memory hierarchies, such as shared memory 132, processor caches, or processor registers.
- the calculation program 134B is a program that implements calculation processing in each processor 133 and data storage processing in the shared memory 132 and storage 134 .
- the control program 134 ⁇ /b>C is a program that controls the calculation server 103 a based on commands sent from the control unit 113 of the management server 101 and sends calculation results of the calculation server 103 a to the management server 101 .
- Such a calculation server 103a executes programs for the processors 133a to 133d to solve combinatorial optimization problems. This program causes the calculation server 103 a to function as the input unit 12 , update unit 14 and output unit 16 .
- the information processing system 100 with such a configuration can be used as a PC cluster.
- a PC cluster is a system that connects a plurality of computers to achieve computational performance that cannot be obtained with a single computer.
- the information processing system 100 can execute parallel calculations by using, for example, MPI (Message Passing Interface), even in a configuration in which memories are distributed among a plurality of calculation servers 103. be.
- MPI Message Passing Interface
- the information processing system 100 may be multiple GPUs connected by high-speed links. In this case, each of the plurality of GPUs executes processing similar to that of the calculation server 103 .
- the present invention is not limited to the above-described embodiments as they are, and can be embodied by modifying the constituent elements without departing from the scope of the present invention at the implementation stage.
- various inventions can be formed by appropriate combinations of the plurality of constituent elements disclosed in the above embodiments. For example, some components may be omitted from all components shown in the embodiments. Furthermore, components across different embodiments may be combined as appropriate.
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
不等式の制約条件の下で0-1組合せ最適化問題を解く。求解装置は、更新部と、出力部とを備える。更新部は、第1変数および第2変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、第1変数および第2変数を交互に更新する。出力部は、終了時刻における複数の要素のそれぞれの第1変数に基づき0-1組合せ最適化問題の解を出力する。単位時間毎の更新処理において、更新部は、1または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第1変数を代入した不等式を満たさない場合、複数の要素のそれぞれの第2変数から、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。
Description
本発明の実施形態は、求解装置、求解方法およびプログラムに関する。
従来、一次不等式による制約の下でイジング問題を解く方法として、スラック変数を利用した定式化をする方法が知られている。この方法は、一次不等式の定数項をWmとした場合、目的関数に、W個(Wは自然数)のスラック変数により表されたペナルティ項を加算して、イジング問題を解く。従って、この方法は、目的関数にN個の決定変数が含まれるならば、N+W個のイジングスピンを含むイジング問題を解く。また、例えば、M個の一次不等式による制限の下でイジング問題を解くのであれば、この方法は、N+{W1+…+Mm+…+WM}個のイジングスピンを含むイジング問題を解く。なお、Mmは、M個の一次不等式のうちのm番目の一次不等式の定数項であり、自然数を表す。
株式のポートフォリオ最適化問題は、リスク回避のために、銘柄の投資割合を決定変数とし、複数の決定変数のそれぞれに重みを乗算した合計を上限値以下とする制約条件の下で、解かれる場合が多い。制約条件は、下記の式のような一次不等式により表される。
xA1・ωA1+xA2・ωA2+…+xSm・ωSm+…≦Wm
xA1・ωA1+xA2・ωA2+…+xSm・ωSm+…≦Wm
この式において、Wmは、上限値を表す。xA1は、銘柄Aの投資割合を例えば1%とする決定変数を表す。xA2は、銘柄Aの投資割合を例えば2%とする決定変数を表す。xSsは、銘柄Sの投資割合をs%とする決定変数を表す。また、ωA1、ωA2およびωSmは、重みを表す。
例えば、株式のポートフォリオ最適化問題は、複数の銘柄を、業種または企業規模毎にグループ化し、グループ毎に上記の制約条件を設定する。例えば、東京証券取引所は、約2000銘柄を取り扱っている。また、東京証券取引所は、これら約2000銘柄を、33種類の業種および3種類の規模に分類している。従って、東京証券所で取り扱われている約2000銘柄は、33業種×3規模の合計99個のグループに分けることができる。
1つの銘柄の投資割合を1%から100%までの100段階で離散化して、東京証券取引所により取り扱われる2000銘柄に対してポートフォリオ最適化問題を解く場合、目的関数に含まれる決定変数の個数は、200000個となる。
これに対して、99個のグループの全てに対して一次不等式による制約条件を設定するとする。上限値(Wm)を1000段階で離散化して、スラッグ変数を利用した定式化をした場合、スラッグ変数の個数は、99000個となる。このため、99個のグループの全てに対して一次不等式による制約条件を設定したポートフォリオ最適化問題を解く場合、イジングスピンの個数が20000+99000=299000個のイジング問題を解かなければならない。
このように、スラック変数を利用した定式化をしてポートフォリオ最適化問題を解く場合、制約条件が無い場合と比較して、例えば東京証券取引所の例では、イジングスピンの数が50%程度増加してしまう。ソルバー等によりイジング問題を解く場合、イジングスピンの数が増加してしまうと、演算量および演算時間が多くなり、求解のためのコストが大きくなってしまう。
H. Goto,"Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network", Sci. Rep. 6, 21686 (2016).
Hayato Goto, Kosuke Tatsumura, Alexander R. Dixon, "Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems",Science Advances, Vol. 5, no. 4, eaav2372, 19 Apr. 2019
土屋和雄,西山岳宏,辻田勝吉, "分岐特性を用いた組合せ最適化問題の近似解法", [online],[2021年5月17日検索],インターネット, <http://www.ynl.t.u-tokyo.ac.jp/project/RobotBrainCREST/publications/pdf/tsuchiya/4_01.pdf>
土屋和雄,西山岳宏,辻田勝吉, "分岐特性を用いた組合せ最適化問題の近似解法 第2報:決定論的アニーリングアルゴリズムの解析", [online],[2021年5月17日検索],インターネット, <http://www.ynl.t.u-tokyo.ac.jp/project/RobotBrainCREST/publications/pdf/tsuchiya/4_02.pdf>
A. Lucas, "Ising formulations of many NP problems", Frontiers in Physics 2, 5 (2014), DOI: 10.3389/fphy.2014.00005
Andrew Lucas, "Ising formulations of many NP problems",Frontiers in Physics, Interdisciplinary Physics, Volume 2, Articles 5, 12 February 2014
本発明が解決しようとする課題は、不等式により表される1または複数の制約条件の下で0-1組合せ最適化問題を少ない演算コストで解くことができる求解装置、求解方法およびプログラムを提供する。
実施形態に係る求解装置は、0-1組合せ最適化問題を、前記0-1組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される1または複数の制約条件の下で解く。前記求解装置は、更新部と、出力部とを備える。前記更新部は、第1変数および第2変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第1変数および前記第2変数を交互に更新する。前記出力部は、前記終了時刻における前記複数の要素のそれぞれの前記第1変数に基づき前記0-1組合せ最適化問題の解を出力する。前記複数の要素は、前記複数の離散変数に対応する。前記第1変数および前記第2変数のそれぞれは、実数により表される。前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第1変数を前記第2変数に基づき更新する。前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第2変数を前記第1変数に基づき更新する。前記単位時間毎の更新処理において、前記更新部は、前記1または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第1変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第2変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。
(0-1組合せ最適化問題)
イジング問題を解くために使われる装置の一例として、イジングマシンが挙げられる。イジングマシンは、イジングモデルの基底状態のエネルギーを計算する。これまで、イジングモデルは、主に強磁性体や相転移現象のモデルとして使われることが多かった。しかし、近年、イジングモデルは、0-1組合せ最適化問題を解くためのモデルとしての利用が増えている。式(1)は、イジングモデルのエネルギーを示す。
イジング問題を解くために使われる装置の一例として、イジングマシンが挙げられる。イジングマシンは、イジングモデルの基底状態のエネルギーを計算する。これまで、イジングモデルは、主に強磁性体や相転移現象のモデルとして使われることが多かった。しかし、近年、イジングモデルは、0-1組合せ最適化問題を解くためのモデルとしての利用が増えている。式(1)は、イジングモデルのエネルギーを示す。
si、sjはスピンを表す。スピンは、+1または-1のいずれかの値をとる2値変数である。siは、i番目のスピンを表す。sjは、j番目のスピンを表す。iおよびjは、1以上、N以下の整数である。Nは、スピンの数を表し、2以上の整数である。hiは、i番目のスピンに作用する局所磁場を表す。Jは、2つのスピン間に作用する力を表す結合係数の行列である。Jは、対角成分が0である実対称行列である。Jijは、Jのi行j列の要素を表す。つまり、Jijは、i番目のスピンと、j番目のスピンとの間に作用する力を表す結合係数である。
イジングマシンは、式(1)により表されるエネルギーEIsingを目的関数とし、エネルギーEIsingを可能な限り小さくする解を算出する。エネルギーEIsingが最小値となるイジングモデルの解(s1、s2、・・・、sN)は、最適解と呼ばれる。ただし、イジングモデルの解は、最適解ではなく、エネルギーEIsingが最小値に近い近似解であってもよい。すなわち、イジング問題は、最適解のみならず、近似解を算出する問題であってもよい。
また、0または1のいずれかの値をとる離散変数(ビット)の2次関数を目的関数とする0-1組合せ最適化問題は、0-1二次計画問題と呼ばれる。離散変数(ビット)は、(1+si)/2の演算を用いることにより、siに変換される。つまり、0-1二次計画問題は、式(1)で表されるイジング問題と等価であるといえる。従って、0-1二次計画問題は、イジング問題に変換し、イジングマシンにより解を算出することが可能である。
イジングマシンは、例えば、量子アニーラ、コヒーレントイジングマシンおよび量子分岐マシン等によりハードウェア実装される。量子アニーラは、超伝導回路を使って量子アニーリングを実現する。コヒーレントイジングマシンは、光パラメトリック発振器で形成されたネットワークの発振現象を利用する。量子分岐マシンは、カー効果を有するパラメトリック発振器のネットワークにおける量子力学的な分岐現象を利用する。これらのハードウェア実装されたイジングマシンは、計算時間の大幅な短縮を実現する可能性がある一方、大規模化および安定的な運用が難しいという課題もある。
イジング問題は、広く普及しているデジタルコンピュータを用いて解を算出することも可能である。デジタルコンピュータは、量子アニーラ、コヒーレントイジングマシンおよび量子分岐マシン等と比べ、大規模化および安定運用が可能である。シミュレーテッドアニーリング(SA)は、デジタルコンピュータでイジング問題を解くためのアルゴリズムの一例である。ただし、シミュレーテッドアニーリングは、それぞれの変数が逐次更新される逐次更新アルゴリズムであるため、並列化による計算処理の高速化は難しい。
非特許文献2には、0-1組合せ最適化問題を解くためのアルゴリズムとして、シミュレーテッド分岐アルゴリズムが提案されている。シミュレーテッド分岐アルゴリズムは、イジングモデルを用いて、デジタルコンピュータによって、規模の大きい0-1組合せ最適化問題を高速に解くことが可能である。シミュレーテッド分岐アルゴリズムは、CPU(Central Processing Unit)、マイクロプロセッサ、GPU(Graphics Processing Unit)、FPGA(Field-Programmable Gate Array)、ASIC(Application Specific Integrated Circuit)、または、これらの組合せの回路等の電子回路によっても、規模の大きい0-1組合せ最適化問題を高速に解くことが可能である。
(シミュレーテッド分岐アルゴリズム)
シミュレーテッド分岐アルゴリズムは、それぞれがN個の要素に対応する変数xiおよび変数yiを用いる。変数xiを第1変数、変数yiを第2変数と呼ぶ場合もある。シミュレーテッド分岐アルゴリズムにおいて、N個の要素のそれぞれは、仮想的な粒子を表す。N個の要素は、イジング問題のN個のスピンに対応する。従って、N個の要素は、組合せ最適化問題のN個の離散変数(ビット)に対応する。変数xiおよび変数yiは、いずれも、実数で表される連続変数である。変数xiは、N個の粒子のうちのi番目の粒子の位置を表す。変数yiは、i番目の粒子の運動量を表す。Nは、2以上の整数である。iは、1以上、N以下の整数を表し、N個の要素のそれぞれを特定するインデックスを表す。
シミュレーテッド分岐アルゴリズムは、それぞれがN個の要素に対応する変数xiおよび変数yiを用いる。変数xiを第1変数、変数yiを第2変数と呼ぶ場合もある。シミュレーテッド分岐アルゴリズムにおいて、N個の要素のそれぞれは、仮想的な粒子を表す。N個の要素は、イジング問題のN個のスピンに対応する。従って、N個の要素は、組合せ最適化問題のN個の離散変数(ビット)に対応する。変数xiおよび変数yiは、いずれも、実数で表される連続変数である。変数xiは、N個の粒子のうちのi番目の粒子の位置を表す。変数yiは、i番目の粒子の運動量を表す。Nは、2以上の整数である。iは、1以上、N以下の整数を表し、N個の要素のそれぞれを特定するインデックスを表す。
係数Dは、予め定められた定数であり、離調(detuning)に相当する。係数p(t)は、ポンピング振幅(pumping amplitude)に相当し、シミュレーテッド分岐アルゴリズムの計算時に更新回数に応じて値が単調増加する。tは、時刻を表す変数である。係数p(t)の初期値は0に設定されていてもよい。係数Kは、予め定められた定数であって、正のカー係数(Kerr coefficient)に相当する。なお、係数Kは、0であってもよい。
式(4)のziは、式(3)の中の小カッコの内の数式を変数xiで偏微分した式である。式(3)の中の小カッコの内の数式は、イジングモデルのエネルギーEIsingに対応する。
cは、係数である。cは、例えば、計算を実行する前に予め定められる定数であってもよい。また、α(t)は、p(t)とともに増加する係数である。
そして、シミュレーテッド分岐アルゴリズムは、p(t)の値を初期値(例えば、0)から所定の値まで増加させた後における変数xiの符号に基づき、スピンsiの値を算出する。シミュレーテッド分岐アルゴリズムは、例えば、xi>0の場合にsgn(xi)=1、xi<0の場合にsgn(xi)=-1となる符号関数を用いて、スピンsiの値を算出する。
(シミュレーテッド分岐アルゴリズムの演算)
シミュレーテッド分岐アルゴリズムは、シンプレクティック・オイラー法を用いて、式(2)、式(3)および式(4)によって与えられる微分方程式を解く。
シミュレーテッド分岐アルゴリズムは、シンプレクティック・オイラー法を用いて、式(2)、式(3)および式(4)によって与えられる微分方程式を解く。
tは、時刻を表す。Δtは、単位時間(時間ステップ、時間刻み幅)を表す。
シミュレーテッド分岐アルゴリズムを実行する場合、デジタルコンピュータまたはFPGA等の電子回路は、式(5)または式(6)のアルゴリズムに基づき、それぞれN個ある変数xiおよび変数yiを初期時刻から単位時間毎に順次に、且つ、変数xiと変数yiとを交互に、更新する。そして、デジタルコンピュータまたはFPGA等の電子回路は、終了時刻におけるN個の変数xiの値を、符号関数を用いて2値化して、N個のスピンの値を出力する。
なお、式(5)および式(6)は、微分方程式との対応関係を示すために、時刻tおよび単位時間Δtを用いて表されている。ただし、シンプレクティック・オイラー法をデジタルコンピュータまたはFPGA等の電子回路で実行する場合、式(5)および式(6)を演算するためのアルゴリズムは、明示的なパラメータとして時刻tおよび単位時間Δtを含まなくてよい。例えば、単位時間Δtを1とする場合、式(5)および式(6)を演算するためのアルゴリズムは、単位時間Δtを含まなくてよい。例えば、明示的なパラメータとして時刻tを含まない場合、式(5)および式(6)を演算するアルゴリズムは、xi(t+Δt)をxi(t)の更新後の値として処理を実行する。すなわち、式(5)および式(6)を演算するアルゴリズムは、“t”を更新前の変数を特定するパラメータ、“t+Δt”を更新後の変数を特定するパラメータとして処理を実行する。以降で説明する式(5)および式(6)を改良したアルゴリズムも同様である。
図1は、シミュレーテッド分岐アルゴリズムにより最適化問題を解いた場合における、変数xiの分岐現象を表す図である。シミュレーテッド分岐アルゴリズムにより最適化問題を解いた場合、系のパラメータが変化することに伴い、安定運動状態が1個のみの系から、安定運動状態が2個の系へと遷移する分岐現象が生じる。図1に示すように、分岐現象が進むと、変数xiは、-1または+1の近傍に集中するが、-1より小さい領域、または、+1より大きい領域にも広がる。
(改良アルゴリズム)
発明者は、上述のシミュレーテッド分岐アルゴリズムを改良して、不等式により表される1または複数の制約条件の下で、0-1組合せ最適化問題を解く改良アルゴリズムを発明した。
発明者は、上述のシミュレーテッド分岐アルゴリズムを改良して、不等式により表される1または複数の制約条件の下で、0-1組合せ最適化問題を解く改良アルゴリズムを発明した。
Gmは、1または複数の制約条件のうちのm番目の制約条件に対応するペナルティ項を表す。1または複数の制約条件のそれぞれは、目的関数に含まれる複数の離散変数を含む不等式により表される。例えば、不等式は、複数の離散変数の一次項を含む一次不等式、または、複数の離散変数の二次項を含む二次不等式である。
ペナルティ項は、対応する不等式を満たす場合に、ハミルトニアンにエネルギーとして0を加え、対応する不等式を満たさない場合に、ハミルトニアンに正のエネルギーを加える。より具体的には、Gmは、制約条件を表す不等式の境界から、N個の要素により特定される位置までの距離が大きい程、大きなエネルギーをハミルトニアンに加える関数である。なお、不等式の境界は、不等式を満たすか満たさないかの境界である。すなわち、不等式の境界は、不等式における不等号を、等号により代えた方程式により表される。
式(8)および式(9)に示す漸化式は、式(5)および式(6)に示す漸化式と比較すると、第2変数(yi(t+Δt))の算出式が異なる。具体的には、式(8)および式(9)に示す漸化式は、第2変数(yi(t+Δt))の算出式に、1または複数の不等式のそれぞれのペナルティ項(Gm)のxiについての偏微分値の合計に、Δtを乗算した値が減算されている点において異なる。
1または複数の不等式のそれぞれの偏微分値は、対応する不等式の境界からN個の要素により特定される位置までの距離における対応する要素の成分に相当する。この偏微分値は、対応する不等式の境界からN個の要素により特定される位置までの距離における対応する要素の成分が大きい程、大きい値となる。
ただし、ペナルティ項は、対応する不等式を満たす場合に、0となる。従って、1または複数の不等式のそれぞれの偏微分値は、N個の要素が不等式を満たしている場合には、0となる。
ここで、シミュレーテッド分岐アルゴリズムは、時刻を表す変数であるtを、単位時間(Δt)毎に増加させながら、複数の第1変数(xi)および複数の第2変数(yi)を交互に更新する。従って、シミュレーテッド分岐アルゴリズムは、単位時間(Δt)毎に、1または複数の制約条件のそれぞれについて、不等式を満たすか否かを判断することができる。
そこで、改良アルゴリズムは、単位時間(Δt)毎に、複数の第2変数(yi)を更新した後に、1または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第1変数(xi)を代入した不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間(Δt)毎に、1または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第1変数を代入した不等式を満たさない場合、複数の第2変数(yi)のそれぞれから、対応する補正値を減算する。つまり、改良アルゴリズムは、単位時間(Δt)毎に、1または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第1変数を代入した不等式を満たさない場合、複数の第2変数(yi)のそれぞれから、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。
これにより、改良アルゴリズムは、複数の第1変数(xi)および複数の第2変数(yi)を単位時間(Δt)毎に、1または複数の制約条件のそれぞれを満たすように、複数の第1変数(xi)および複数の第2変数(yi)を補正することができる。この結果、改良アルゴリズムは、終了時刻(t=T)における複数の第1変数(xi)が1または複数の制約条件の全てを満たす確率を高くすることができる。
このような改良アルゴリズムは、不等式により表される1または複数の制約条件の下で、0-1組合せ最適化問題を解く場合であっても、制約条件を設定しない場合と同一の変数の数により、解を得ることができる。従って、このような改良アルゴリズムは、スラッグ変数を設けずに、簡易に、不等式により表される1または複数の制約条件の下で0-1組合せ最適化問題を解くことができる。
(制約条件が一次不等式であるのペナルティ項および偏微分値の第1例)
1または複数の制約条件の一つは、一次不等式により表される。例えば、改良アルゴリズムは、M個(Mは1以上の整数)の一次不等式の制約の下で、0-1組合せ最適化問題を解く。この場合において、1または複数の制約条件のうちのm番目の制約条件は、式(10)により表される。
1または複数の制約条件の一つは、一次不等式により表される。例えば、改良アルゴリズムは、M個(Mは1以上の整数)の一次不等式の制約の下で、0-1組合せ最適化問題を解く。この場合において、1または複数の制約条件のうちのm番目の制約条件は、式(10)により表される。
式(10)において、mは、1以上、M以下の整数である。式(10)において、xiは、i番目の離散変数である。式(10)において、xjは、j番目の離散変数である。ωm,iは、i番目の離散変数に乗算される係数である。Wmは、予め定められた定数である。
Amは、m番目の制約条件に対して予め定められた係数である。kは、1以上の整数である。例えば、kは、2であってもよい。
Smは、一次不等式における定数項(Wm)を除く関数における複数の離散変数のそれぞれに、対応する第1変数(xi)を代入して得られる値である。具体的には、Smは、式(8)の漸化式を適用する場合には、式(12)のように表される。また、Smは、式(9)の漸化式を適用する場合には、式(13)のように表される。
そこで、制約条件が一次不等式で表される場合、改良アルゴリズムは、単位時間(Δt)毎に、複数の第2変数(yi)を更新した後に、複数の離散変数のそれぞれに対応する第1変数(xi)を代入した一次不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間(Δt)毎に、対応する一次不等式を満たさない場合、複数の第2変数(yi)のそれぞれから、式(14)に表される偏微分値に単位時間を表すΔtを乗算した補正値を、減算する。
このような処理を実行する改良アルゴリズムは、制約条件を設定しない場合と同一の変数の数により、一次不等式により表される制約条件の下で0-1組合せ最適化問題を簡易に解くことができる。
(制約条件が一次不等式であるのペナルティ項および偏微分値の第2例)
式(14)に示す偏微分値は、Amが大きい場合、計算誤差が蓄積して不安定となる可能性がある。そこで、改良アルゴリズムは、制約条件が一次不等式である場合、式(15)により表される偏微分値を用いてもよい。すなわち、改良アルゴリズムは、単位時間(Δt)毎に、対応する一次不等式を満たさない場合、複数の第2変数(yi)のそれぞれから、式(15)に表される偏微分値に単位時間を表すΔtを乗算した補正値を、減算してもよい。
式(14)に示す偏微分値は、Amが大きい場合、計算誤差が蓄積して不安定となる可能性がある。そこで、改良アルゴリズムは、制約条件が一次不等式である場合、式(15)により表される偏微分値を用いてもよい。すなわち、改良アルゴリズムは、単位時間(Δt)毎に、対応する一次不等式を満たさない場合、複数の第2変数(yi)のそれぞれから、式(15)に表される偏微分値に単位時間を表すΔtを乗算した補正値を、減算してもよい。
Bmは、0以上1以下の値であり、m番目の制約条件に対して予め定められた係数である。Emは、一次不等式における定数項(Wm)を除く関数における複数の離散変数のそれぞれに、対応する第2変数(yi)を代入して得られる値である。具体的には、Emは、式(8)または式(9)の漸化式を適用する場合には、式(16)のように表される。
Emは、複数の要素により特定される方向の、一次不等式の境界に対する傾きを表す。従って、式(15)の偏微分値を用いた場合、補正値は、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分と、複数の要素により特定される方向の、一次不等式の境界に対する傾きの対応する要素の成分とを加算した値に応じた値となる。これにより、改良アルゴリズムは、Amが大きいために距離に応じた成分が不安定となる場合であっても、補正値を安定化させることができる。
式(17)において、xiは、i番目の離散変数である。式(17)において、xjは、j番目の離散変数である。Qi,jは、N行、N列の半正定値行列に含まれるi行、j列の値である。qは、予め定められた定数である。
Aは、二次不等式による制約条件に対して予め定められた係数である。kは、1以上の整数である。例えば、kは、2であってもよい。
Sは、二次不等式における定数項(q)を除く関数における複数の離散変数のそれぞれに、対応する第1変数(xi,xj)を代入して得られる値である。具体的には、Sは、式(8)の漸化式を適用する場合には、式(19)のように表される。また、Sは、式(9)の漸化式を適用する場合には、式(20)のように表される。
そこで、制約条件が二次不等式で表される場合、改良アルゴリズムは、単位時間(Δt)毎に、複数の第2変数(yi)を更新した後に、複数の離散変数のそれぞれに対応する第1変数(xi,xj)を代入した二次不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間(Δt)毎に、対応する二次不等式を満たさない場合、複数の第2変数(yi)のそれぞれから、式(21)または式(22)に表される偏微分値に単位時間を表すΔtを乗算した補正値を、減算する。
イジングモデルは、二次不等式による制約条件の下で0-1組合せ最適化問題を解くことは困難である。これに対して、改良アルゴリズムは、二次不等式により表される制約条件の下で、0-1組合せ最適化問題を、簡易に解を得ることができる。
なお、改良アルゴリズムは、1または複数の制約条件の中に、一次不等式による制約条件と、二次不等式による制約条件との両者が含まれていても、0-1組合せ最適化問題を解くことができる。また、改良アルゴリズムは、1または複数の制約条件に一次不等式による制約条件のみが含まれていてもよい。また、改良アルゴリズムは、1つの二次不等式による制約条件のみの下で、0-1組合せ最適化問題を解いてもよい。
(機能ブロック構成)
図2は、本実施形態に係る求解装置10の機能構成を示す図である。
図2は、本実施形態に係る求解装置10の機能構成を示す図である。
求解装置10は、改良アルゴリズムを用いて、0-1組合せ最適化問題を不等式により表される1または複数の制約条件の下で解く。求解装置10は、コンピュータ等の情報処理装置、ネットワークを介して複数のコンピュータまたはサーバが相互に通信をして構成されるコンピュータシステム、または、複数台のコンピュータが連携して情報処理を実行するPCクラスタ等により実現される。また、求解装置10は、CPU、マイクロプロセッサ、GPU、FPGAまたはASIC、または、これらの組合せの回路等の電子回路によって実現される。
求解装置10は、機能構成として、入力部12と、更新部14と、出力部16とを備える。
入力部12は、0-1組合せ最適化問題の目的関数を定義するための情報(例えば、N、J、h)、および、改良アルゴリズムを実行するために必要な係数を表す情報(例えば、D、c、Δt、T,p(t)、α(t))を外部装置から受け取る。さらに、入力部12は、1または複数の制約条件を定義するための情報(例えば、ωm,i,Wm,Qi,j,q)、および、制約処理を実行するために必要な係数(例えば、Am,k,Bm,A)を表す情報を外部装置から受け取る。そして、入力部12は、受け取ったこれらの情報を更新部14に与える。
更新部14は、改良アルゴリズムを用いて、第1変数(xi)および第2変数(yi)が対応付けられた複数の要素のそれぞれについて、初期時刻(t=0)から終了時刻(t=T)まで単位時間(Δt)毎に順次に、第1変数(xi)および第2変数(yi)を交互に更新する。
出力部16は、終了時刻(t=T)における複数の要素のそれぞれの第1変数(xi)に基づき、0-1組合せ最適化問題の解を出力する。例えば、出力部16は、終了時刻における複数の要素のそれぞれについて、第1変数(xi)を予め設定されたしきい値により2値化した離散変数の値を算出する。そして、出力部16は、算出した複数の離散変数の値を0-1組合せ最適化問題の解として出力する。
ここで、複数の要素は、0-1組合せ最適化問題の複数の離散変数に対応する。また、第1変数(xi)および第2変数(yi)のそれぞれは、実数により表される。
そして、単位時間毎の更新処理において、更新部14は、複数の要素のそれぞれについて、第1変数(xi)を第2変数(yi)に基づき更新する。また、単位時間毎の更新処理において、更新部14は、複数の要素のそれぞれについて、第2変数(yi)を第1変数(xi)に基づき更新する。
例えば、単位時間毎の更新処理において、更新部14は、複数の要素のそれぞれについて、第1変数(xi)を更新した後に第2変数(yi)を更新する。これに代えて、単位時間毎の更新処理において、更新部14は、複数の要素のそれぞれについて、第2変数(yi)を更新した後に第1変数(xi)を更新してもよい。
さらに、単位時間毎の更新処理において、更新部14は、第2変数(yi)を更新した後に、0-1組合せ最適化問題の解が1または複数の制約条件を満たすように、制約処理を実行する。更新部14は、1または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第1変数(xi)を代入した不等式を満たさない場合、複数の要素のそれぞれの第2変数(yi)から、対応する不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値(rm,i)を減算する。
単位時間毎の更新処理において、1または複数の制約条件のそれぞれについて制約処理を実行することにより、更新部14は、複数の第1変数(xi)が制約条件を満たしていない場合、複数の第1変数(xi)が制約条件を満たす方向に変化するように、複数の要素に対応する複数の第2変数(yi)を補正することができる。これにより、更新部14は、終了時刻(t=T)における複数の第1変数(xi)が1または複数の制約条件の全てを満たす確率を高くすることができる。
なお、単位時間毎の更新処理において、更新部14は、第1変数(xi)を更新した後に第2変数(yi)を更新する場合には、複数の要素のそれぞれについて、第2変数(yi)を更新した後に第2変数(yi)から補正値(rm,i)を減算する。また、単位時間毎の更新処理において、更新部14は、第2変数(yi)を更新した後に第1変数(xi)を更新する場合には、複数の要素のそれぞれについて、第2変数(yi)を更新した後、且つ、第1変数(xi)を更新する前に、第2変数(yi)から補正値(rm,i)を減算する。
(処理フロー)
図3は、更新部14の処理の流れの第1例を示すフローチャートである。更新部14は、例えば、図3に示す流れで処理を実行する。
図3は、更新部14の処理の流れの第1例を示すフローチャートである。更新部14は、例えば、図3に示す流れで処理を実行する。
まず、S101において、更新部14は、0-1組合せ最適化問題を解くためのパラメータを設定する。具体的には、更新部14は、N×N個の結合係数を含む行列であるJ、および、N個の局所磁場を表す局所磁場係数を含む配列であるhを設定する。さらに、更新部14は、係数であるD、係数であるc、単位時間を表すΔt、終了時刻を表すT、関数であるp(t)、および、関数であるα(t)を設定する。p(t)およびα(t)は、t=初期時刻(例えば0)で0、t=終了時刻(T)で1となる増加関数である。更新部14は、J、hを入力部12からの受け取った情報に応じて設定する。更新部14は、D、c、Δt、T、p(t)およびα(t)を、入力部12から受け取った値に応じて設定してもよいし、予め決定されており変更できない値を設定してもよい。
続いて、S102において、更新部14は、制約条件に関する情報を設定する。具体的には、更新部14は、1または複数の制約条件を定義するための情報、および、制約処理を実行するために必要な係数を表す情報を設定する。例えば、更新部14は、制約条件が一次不等式により表される場合、変数に乗算される係数であるωm,i、定数項であるWm、制約処理を実行するために必要な係数であるAm(またはAmおよびBm)、並びに、kを設定する。また、例えば、制約条件が二次不等式により表される場合、更新部14は、変数に乗算される係数であるQi,jおよび定数項であるq、制約処理を実行するために必要な係数であるAおよびkを設定する。なお、更新部14は、制約処理を実行するために必要な係数であるAm、Bm、kおよびAを、入力部12から受け取ったパラメータに応じて設定してもよいし、予め決定されており変更できない値を設定してもよい。
続いて、S103において、更新部14は、変数を初期化する。具体的には、更新部14は、時刻を表す変数であるtを初期時刻(例えば、0)に初期化する。さらに、更新部14は、N個の第1変数(x1(t)~xN(t))のそれぞれおよびN個の第2変数(y1(t)~yN(t))のそれぞれに、ユーザから受け取った初期値、予め定められた固定値、または、乱数を代入する。
続いて、更新部14は、S104とS116との間のループ処理を、tがTより大きくなるまで繰り返す。1回のループ処理において、更新部14は、対象時刻(t+Δt)におけるN個の第1変数(x1(t+Δt)~xN(t+Δt))を、直前時刻(t)におけるN個の第1変数(x1(t)~xN(t))、および、直前時刻(t)におけるN個の第2変数(y1(t)~yN(t))に基づき算出する。また、1回のループ処理において、更新部14は、対象時刻(t+Δt)におけるN個の第2変数(y1(t+Δt)~yN(t+Δt))を、対象時刻(t+Δt)におけるN個の第1変数(x1(t+Δt)~xN(t+Δt))および直前時刻(t)におけるN個の第2変数(y1(t)~yN(t))に基づき算出する。
なお、直前時刻(t)は、対象時刻(t+Δt)より単位時間(Δt)前の時刻である。すなわち、更新部14は、S104とS116との間のループ処理を繰り返すことにより、N個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を、初期時刻(t=0)から終了時刻(t=T)まで単位時間(Δt)毎に順次に更新する。
続いて、更新部14は、S105とS107との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。iは、1からNまでの整数であり、N個の要素のうちの処理対象を表すインデックスである。N個の要素のそれぞれは、第1変数(xi(t))および第2変数(yi(t))が対応付けられる。S105とS107との間のループ処理において、更新部14は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S106において、更新部14は、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、対象要素の直前時刻(t)における第1変数(xi(t))に、対象要素の直前時刻(t)における第2変数(yi(t))と予め定められた定数(D)と単位時間(Δt)とを乗算した値を加算することにより算出する。具体的には、更新部14は、式(23)を算出する。
すなわち、更新部14は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、対象要素の直前時刻(t)における第1変数(xi(t))と、対象要素の直前時刻(t)における第2変数(yi(t))とに基づき更新する。
更新部14は、S105とS107との間のループ処理をN回実行した場合、処理をS108に進める。
続いて、更新部14は、S108とS113との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。
S109において、更新部14は、N個の要素のそれぞれの対象時刻(t+Δt)における第1変数(x1(t+Δt)~xN(t+Δt))と、対象要素とN個の要素のそれぞれとの組毎に0-1組合せ最適化問題により予め定められる作用係数と、に基づき更新値(zi(t+Δt))を算出する。作用係数は、Jに含まれる結合係数およびhに含まれる局所磁場係数である。具体的には、更新部14は、式(24)を算出する。
続いて、S111において、更新部14は、時間経過に従って増加する関数であるp(t+Δt)に基づき定まる値に、対象要素の対象時刻(t+Δt)における第1変数(x1(t+Δt))を乗算した時間発展値(gi(t+Δt))を算出する。具体的には、更新部14は、式(26)を算出する。
続いて、S112において、更新部14は、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δ))を、対象要素の直前時刻(t)における第2変数(yi(t))に、時間発展値(gi(t+Δt))と外力(fi(t+Δt))とを加算した値に単位時間(Δt)を乗算した値を、加算することにより、算出する。具体的には、更新部14は、式(27)を算出する。
更新部14は、以上のようなS108とS113との間のループ処理をN回実行することにより、N個の要素のそれぞれについて、対象時刻(t+Δt)における第2変数(yi(t+Δt))を、対象時刻(t+Δt)におけるN個の第1変数(x1(t+Δt)~xN(t+Δt))と、対象要素の直前時刻(t)におけるに第2変数(yi(t))とに基づき更新する。
更新部14は、S108とS113との間のループ処理をN回実行した場合、処理をS114に進める。
S114において、更新部14は、1または複数の制約条件に基づく制約処理を実行する。なお、制約処理については、図5~図7のフローチャートを用いて詳細を後述する。S114の処理を終えると、更新部14は、処理をS115に進める。
S115において、更新部14は、直前時刻(t)および対象時刻(t+Δt)のそれぞれに単位時間(Δt)を加算して、直前時刻(t)および対象時刻(t+Δt)を更新する。S116において、更新部14は、S105からS115までの処理を、tが終了時刻(T)を超えるまで繰り返す。そして、更新部14は、tが終了時刻(T)より大きくなった場合、本フローを終了する。
そして、出力部16は、N個の要素のそれぞれについて、終了時刻(t=T)における第1変数(xi(T))の符号に応じて、対応するスピンの値を算出する。例えば、出力部16は、終了時刻(t=T)における第1変数(xi(T))の符号が負である場合、対応するスピンを-1とし、正である場合、対応するスピンを+1とする。そして、出力部16は、算出した複数のスピンの値、または、算出した複数のスピンの値を離散変数に変換した値を組合せ最適化問題の解として出力する。
以上のS101~S116の処理を実行することにより、更新部14は、改良アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
図4は、更新部14の処理の流れの第2例を示すフローチャートである。更新部14は、改良アルゴリズムを用いて0-1組合せ最適化問題を解く場合、図3に示す流れに代えて、図4に示す流れで処理を実行してもよい。
まず、S201、S202およびS203において、更新部14は、図3に示す第1例のS101、S102およびS103と同一の処理を実行する。
続いて、更新部14は、S204とS216との間のループ処理を、tがTより大きくなるまで繰り返す。1回のループ処理において、更新部14は、対象時刻(t+Δt)におけるN個の第2変数(y1(t+Δt)~yN(t+Δt))を、直前時刻(t)におけるN個の第1変数(x1(t)~xN(t))および直前時刻(t)におけるN個の第2変数(y1(t)~yN(t))に基づき算出する。また、1回のループ処理において、更新部14は、対象時刻(t+Δt)におけるN個の第1変数(x1(t+Δt)~xN(t+Δt))を、直前時刻(t)におけるN個の第1変数(x1(t)~xN(t))、および、対象時刻(t+Δt)におけるN個の第2変数(y1(t+Δt)~yN(t+Δt))に基づき算出する。
続いて、更新部14は、S205とS210との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。S205とS210との間のループ処理において、更新部14は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S206において、更新部14は、N個の要素のそれぞれの直前時刻(t)における第1変数(x1(t)~xN(t))と、対象要素とN個の要素のそれぞれとの組毎に0-1組合せ最適化問題により予め定められる作用係数と、に基づき更新値(zi(t))を算出する。具体的には、更新部14は、式(28)を算出する。
続いて、S208において、更新部14は、時間経過に従って増加する関数であるp(t)に基づき定まる値に、対象要素の直前時刻(t)における第1変数(x1(t))を乗算した時間発展値(gi(t))を算出する。具体的には、更新部14は、式(30)を算出する。
続いて、S209において、更新部14は、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δ))を、対象要素の直前時刻(t)における第2変数(yi(t))に、時間発展値(gi(t))と外力(fi(t))とを加算した値に単位時間(Δt)を乗算した値を、加算することにより、算出する。具体的には、更新部14は、式(31)を算出する。
更新部14は、以上のようなS205とS210との間のループ処理をN回実行することにより、N個の要素のそれぞれについて、対象時刻(t+Δt)における第2変数(yi(t+Δt))を、直前時刻(t)におけるN個の第1変数(x1(t)~xN(t))と、対象要素の直前時刻(t)におけるに第2変数(yi(t))とに基づき更新する。
更新部14は、S205とS210との間のループ処理をN回実行した場合、処理をS211に進める。
S211において、更新部14は、1または複数の制約条件に基づく制約処理を実行する。なお、制約処理については、図5~図7のフローチャートを用いて詳細を後述する。S211の処理を終えると、更新部14は、処理をS212に進める。
続いて、更新部14は、S212とS214との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。
S213において、更新部14は、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、対象要素の直前時刻(t)における第1変数(xi(t))に、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δt))と予め定められた定数(D)と単位時間(Δt)とを乗算した値を加算することにより算出する。具体的には、更新部14は、式(32)を算出する。
すなわち、更新部14は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、対象要素の直前時刻(t)における第1変数(xi(t))と、対象要素の直前時刻(t)における第2変数(yi(t))とに基づき更新する。
更新部14は、S212とS214との間のループ処理をN回実行した場合、処理をS215に進める。
S215において、更新部14は、直前時刻(t)および対象時刻(t+Δt)のそれぞれに単位時間(Δt)を加算して、直前時刻(t)および対象時刻(t+Δt)を更新する。S216において、更新部14は、S212からS215までの処理を、tが終了時刻(T)を超えるまで繰り返す。そして、更新部14は、tが終了時刻(T)より大きくなった場合、本フローを終了する。
そして、出力部16は、N個の要素のそれぞれについて、終了時刻(t=T)における第1変数(xi(T))の符号に応じて、対応するスピンの値を算出する。例えば、出力部16は、終了時刻(t=T)における第1変数(xi(T))の符号が負である場合、対応するスピンを-1とし、正である場合、対応するスピンを+1とする。そして、出力部16は、算出した複数のスピンの値、または、算出した複数のスピンの値を離散変数に変換した値を組合せ最適化問題の解として出力する。
以上のS201~S216の処理を実行することにより、更新部14は、改良アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
(制約処理の第1例)
図5は、制約処理の第1例の処理の流れを示すフローチャートである。制約条件として、1または複数の制約条件としてM個の一次不等式が設定されている場合、更新部14は、図3のS114および図4のS211において、図5に示す流れで処理を実行する。
図5は、制約処理の第1例の処理の流れを示すフローチャートである。制約条件として、1または複数の制約条件としてM個の一次不等式が設定されている場合、更新部14は、図3のS114および図4のS211において、図5に示す流れで処理を実行する。
まず、更新部14は、S301とS308との間のループ処理を、m=1からi=Mまでmを1ずつインクリメントしながら繰り返す。なお、mは、M個の一次不等式のうちの処理対象となる一次不等式を表すインデックスである。
S301とS308とのループ内において、まず、S302において、更新部14は、制約値(Sm)を算出する。具体的には、図3のS114の制約処理において、更新部14は、式(33)の演算を実行する。また、図4のS211の制約処理において、更新部14は、式(34)の演算を実行する。
なお、制御値であるSmは、対応する一次不等式における定数項(Wm)を除く関数に含まれる複数の離散変数のそれぞれに、対応する第1変数(xi)を代入し得られる値である。
続いて、S303において、更新部14は、制約値(Sm)が対応する一次不等式における定数項(Wm)以下であるか否かを判断する。すなわち、更新部14は、複数の第1変数(xi)が対応する一次不等式を満たすか否かを判断する。
更新部14は、制約値(Sm)が対応する一次不等式における定数項(Wm)以下である場合、すなわち、複数の第1変数(xi)が対応する一次不等式を満たす場合(S303のYes)、処理をS308に進める。更新部14は、制約値(Sm)が対応する一次不等式における定数項(Wm)以下ではない場合、すなわち、複数の第1変数(xi)が対応する一次不等式を満たさない場合(S303のNo)、処理をS304に進める。
続いて、更新部14は、S304とS307との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。S304とS307との間のループ処理において、更新部14は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S304とS307とのループ内において、まず、S305において、更新部14は、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値(rm,i)を算出する。具体的には、更新部14は、式(35)の演算を実行する。
更新部14は、S304とS307との間のループ処理をN回実行した場合、処理をS308に進める。
そして、更新部14は、S301とS308との間のループ処理をM回実行した場合、すなわち、M個の一次不等式の全てについて処理を実行した場合、本フローを終了する。
以上の処理を実行することにより、更新部14は、M個の一次不等式のそれぞれについて、複数の第1変数(xi)が対応する一次不等式を満たさない場合、複数の要素のそれぞれの第2変数(yi)から、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値(rm,i)を減算することができる。これにより、更新部14は、複数の第1変数(xi)が対応する一次不等式を満たしていない場合、複数の第1変数(xi)が対応する一次不等式を満たす方向に変化するように、複数の第2変数(yi)を補正することができる。この結果、更新部14は、終了時刻(t=T)における複数の第1変数(xi)が1または複数の制約条件の全てを満たす確率を高くすることができる。
(制約処理の第2例)
図6は、制約処理の第2例の処理の流れを示すフローチャートである。制約条件として、1または複数の制約条件としてM個の一次不等式が設定されている場合、更新部14は、図3のS114および図4のS211において、図6に示す流れで処理を実行してもよい。
図6は、制約処理の第2例の処理の流れを示すフローチャートである。制約条件として、1または複数の制約条件としてM個の一次不等式が設定されている場合、更新部14は、図3のS114および図4のS211において、図6に示す流れで処理を実行してもよい。
まず、更新部14は、S401とS409との間のループ処理を、m=1からi=Mまでmを1ずつインクリメントしながら繰り返す。
S401とS409とのループ内において、まず、S402において、更新部14は、制約値(Sm)を算出する。S402の処理は、図5のS302の処理と同一である。
続いて、S403において、更新部14は、制約値(Sm)が対応する一次不等式における定数項(Wm)以下であるか否かを判断する。S403の処理は、図5のS303の処理と同一である。更新部14は、制約値(Sm)が対応する一次不等式における定数項(Wm)以下である場合(S403のYes)、処理をS409に進める。更新部14は、制約値(Sm)が対応する一次不等式における定数項(Wm)以下ではない場合(S403のNo)、処理をS404に進める。
続いて、更新部14は、S405とS408との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。S405とS408との間のループ処理において、更新部14は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S405とS408とのループ内において、まず、S406において、更新部14は、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分、および、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた補正値(rm,i)を算出する。具体的には、更新部14は、式(38)の演算を実行する。
続いて、S407において、更新部14は、対象要素の第2変数(yi)から、補正値(rm,i)を減算する。S407の処理は、図5のS306の処理と同一である。
更新部14は、S405とS408との間のループ処理をN回実行した場合、処理をS409に進める。
そして、更新部14は、S401とS409との間のループ処理をM回実行した場合、すなわち、M個の一次不等式の全てについて処理を実行した場合、本フローを終了する。
以上の処理を実行することにより、更新部14は、M個の一次不等式のそれぞれについて、複数の第1変数(xi)が対応する一次不等式を満たさない場合、複数の要素のそれぞれの第2変数(yi)から、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分、および、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた補正値(rm,i)を減算することができる。特に、更新部14は、補正値(rm,i)に、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた値を加えることにより、単位時間(Δt)毎の更新処理において、複数の第1変数(xi)および複数の第2変数(yi)が不安定に変動することを抑制することができる。
これにより、更新部14は、複数の第1変数(xi)が対応する一次不等式を満たしていない場合、複数の第1変数(xi)が対応する一次不等式を満たす方向に安定して変化するように、複数の要素に対応する複数の第2変数(yi)を補正することができる。この結果、更新部14は、終了時刻(t=T)における複数の第1変数(xi)が1または複数の制約条件の全てを満たす確率を高くすることができる。
(制約処理の第3例)
図7は、制約処理の第3例の処理の流れを示すフローチャートである。1または複数の制約条件のうちの一つに二次不等式が設定されている場合、更新部14は、二次不等式については、図5のS302~S307または図6のS402~S408の処理に代えて、図7に示すS501~S506の処理を実行する。
図7は、制約処理の第3例の処理の流れを示すフローチャートである。1または複数の制約条件のうちの一つに二次不等式が設定されている場合、更新部14は、二次不等式については、図5のS302~S307または図6のS402~S408の処理に代えて、図7に示すS501~S506の処理を実行する。
まず、S501において、更新部14は、制約値(S)を算出する。具体的には、図3のS114の制約処理において、更新部14は、式(39)の演算を実行する。また、図4のS211の制約処理において、更新部14は、式(40)の演算を実行する。
なお、制御値であるSは、二次不等式における定数項(q)を除く関数に含まれる複数の離散変数のそれぞれに、対応する第1変数(xi)を代入して得られる値である。
続いて、S502において、更新部14は、制約値(S)が二次不等式における定数項(q)以下であるか否かを判断する。すなわち、更新部14は、複数の第1変数(xi)が二次不等式を満たすか否かを判断する。
更新部14は、制約値(S)が二次不等式における定数項(q)以下である場合、すなわち、複数の第1変数(xi)が対応する二次不等式を満たす場合(S502のYes)、本フローを終了する。更新部14は、制約値(S)が二次不等式における定数項(q)以下ではない場合、すなわち、複数の第1変数(xi)が対応する二次不等式を満たさない場合(S502のNo)、処理をS503に進める。
続いて、更新部14は、S503とS506との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。S503とS506との間のループ処理において、更新部14は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S503とS506とのループ内において、まず、S504において、更新部14は、二次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値(ri)を算出する。具体的には、図3のS114の制約処理において、更新部14は、式(41)の演算を実行する。また、図4のS211の制約処理において、更新部14は、式(42)の演算を実行する。
更新部14は、S503とS506との間のループ処理をN回実行した場合、本フローを終了する。
以上の処理を実行することにより、更新部14は、二次不等式について、複数の第1変数(xi)が対応する二次不等式を満たさない場合、複数の要素のそれぞれの第2変数(yi)から、二次不等式の境界から複数の要素により特定される位置までの距離の対象要素の成分に応じた補正値(ri)を減算することができる。これにより、更新部14は、複数の第1変数(xi)が対応する二次不等式を満たしていない場合、複数の第1変数(xi)が対応する二次不等式を満たす方向に変化するように、複数の第2変数(yi)を補正することができる。この結果、更新部14は、終了時刻(t=T)における複数の第1変数(xi)が1または複数の制約条件の全てを満たす確率を高くすることができる。
以上のように、本実施形態に係る求解装置10は、シミュレーテッド分岐アルゴリズムにより0-1組合せ最適化問題を解く。本実施形態に係る求解装置10は、単位時間(Δt)毎の更新処理において、1または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第1変数を代入した不等式を満たさない場合、複数の要素のそれぞれの第2変数(yi)から、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値(ri)を減算する。
これにより、本実施形態に係る求解装置10は、不等式により表される1または複数の制約条件の下で、0-1組合せ最適化問題を解く場合であっても、制約条件を設定しない場合と同一の変数の数により、解を得ることができる。従って、本実施形態に係る求解装置10によれば、スラッグ変数を設けずに、簡易に、不等式により表される1または複数の制約条件の下で0-1組合せ最適化問題を解くことができる。この結果、本実施形態に係る求解装置10は、変数の途中経過等を記憶するメモリ量を少なくし、計算量および計算時間を少なくすることができる。
求解装置10は、一回の列番号および行番号の指定により、ωm,iおよびWmを同時に読み出すことができ、容易に処理を実行することができる。
(システム構成)
図8は、情報処理システム100の構成を示す図である。改良アルゴリズムは、例えば図8に示すような情報処理システム100により実行させることが可能である。これにより、情報処理システム100は、大規模な組合せ最適化問題を、並列処理により高速に解くことができる。
図8は、情報処理システム100の構成を示す図である。改良アルゴリズムは、例えば図8に示すような情報処理システム100により実行させることが可能である。これにより、情報処理システム100は、大規模な組合せ最適化問題を、並列処理により高速に解くことができる。
情報処理システム100は、管理サーバ101と、ネットワーク102と、複数の計算サーバ103(103a~103c)(情報処理装置)と、複数のケーブル104(104a~104c)と、スイッチ105を備える。また、端末装置106は、情報処理システム100と通信可能である。管理サーバ101、複数の計算サーバ103(103a~103c)、端末装置106は、ネットワーク102を介して互いにデータ通信をする。ネットワーク102は、例えば、複数のコンピュータネットワークが相互に接続されたインターネットである。ネットワーク102は、通信媒体として有線、無線、または、これらの組合せであってよい。
また、複数の計算サーバ103(103a~103c)は、それぞれケーブル104(104a~104c)を介してスイッチ105に接続される。複数のケーブル104(104a~104c)およびスイッチ105は、計算サーバ103と計算サーバ103と間のインターコネクトを形成する。複数の計算サーバ103(103a~103c)のそれぞれは、インターコネクトを介して互いにデータ通信をすることも可能である。スイッチ105は、例えば、Infinibandのスイッチでる。ケーブル104a~104cは、例えば、Infinibandのケーブルである。ただし、スイッチ105およびケーブル104は、有線LANのスイッチ/ケーブルであってもよい。ケーブル104およびスイッチ105で使われる通信規格および通信プロトコルについては、特に問わない。端末装置106は、例えば、ノートPC、デスクトップPC、スマートフォン、タブレットまたは車載端末装置である。
改良アルゴリズムを用いた組合せ最適化問題の求解処理は、並列的および分散的に実行可能である。従って、複数の計算サーバ103(103a~103c)のそれぞれおよび/または計算サーバ103(103a~103c)のプロセッサは、改良アルゴリズムを用いた組合せ最適化問題の求解処理における一部の計算処理の一部のステップを分担して実行してもよいし、異なる変数に対する同一の計算処理を並列的に実行してもよい。この場合、管理サーバ101は、例えば、ユーザによって入力された組合せ最適化問題を各計算サーバ103に処理可能な形式に変換し、各計算サーバ103に実行させる。そして、管理サーバ101は、各計算サーバ103から計算結果を取得し、集約した計算結果を組合せ最適化問題の解に変換する。
図8には、3台の計算サーバ103(103a~103c)が示されている。しかし、情報処理システム100に含まれる計算サーバ103の台数は、3台に限定されない。例えば、情報処理システム100は、1台であっても、2台であっても、4台以上であってもよい。また、情報処理システム100は、含んでいる複数の計算サーバ103のうちの一部を用いて、組合せ最適化問題の求解の処理を実行してもよい。計算サーバ103は、どのような種類の情報処理装置であってもよい。例えば、計算サーバ103は、データセンターに設置されたサーバであってもよいし、オフィスに設置されたデスクトップPCであってもよい。また、計算サーバ103は、異なるローケーションに設置された複数の種類のコンピュータであってもよい。例えば、計算サーバ103は、汎用的なコンピュータであってもよいし、専用の電子回路または、これらの組合せであってもよい。
図9は、管理サーバ101の構成を示す図である。管理サーバ101は、例えば、中央演算処理装置(CPU)とメモリとを含むコンピュータである。管理サーバ101は、プロセッサ110と、記憶部114と、通信回路115と、入力回路116と、出力回路117とを備える。プロセッサ110、記憶部114、通信回路115、入力回路116、出力回路117は、互いにバス120を介して接続される。プロセッサ110は、内部の機能構成として、管理部111と、変換部112と、制御部113を含む。
プロセッサ110は、演算を実行し、管理サーバ101の制御を行う電子回路である。プロセッサ110は、例えば、CPU、マイクロプロセッサ、ASIC、FPGA、PLDまたはこれらの組合せであってよい。管理部111は、ユーザの端末装置106を介して管理サーバ101の操作を行うためのインタフェースを提供する。管理部111が提供するインタフェースは、例えば、API、CLIまたはウェブページ等である。例えば、ユーザは、管理部111を介して組合せ最適化問題の情報の入力を行ったり、計算された組合せ最適化問題の解の閲覧および/またはダウンロードを行ったりする。変換部112は、組合せ最適化問題に関するパラメータを入力して、入力したパラメータを各計算サーバ103が処理可能な形式に変換する。制御部113は、各計算サーバ103に制御指令を送信する。制御部113が各計算サーバ103から計算結果を取得した後、変換部112は、複数の計算結果を集約し、組合せ最適化問題の解に変換し、組合せ最適化問題の解を出力する。
記憶部114は、管理サーバ101のプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。プログラムは、OSとアプリケーションの両方を含む。記憶部114は、揮発性メモリ、不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリは、例えば、DRAMまたはSRAM等である。不揮発性メモリは、NANDフラッシュメモリ、NORフラッシュメモリ、ReRAMまたはMRAM等である。また、記憶部114は、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置等であってもよい。
通信回路115は、ネットワーク102に接続された各装置との間でデータの送受信を行う。通信回路115は、例えば、有線LANのNIC(Network Interface Card)である。ただし、通信回路115は、無線LANなど、その他の種類の通信回路であってもよい。入力回路116は、管理サーバ101へのデータ入力を実現する。入力回路116は、外部ポートとして、例えば、USB、PCI-Expressなどを含む。操作装置118は、入力回路116に接続される。操作装置118は、ユーザが管理サーバ101に情報を入力するための装置である。操作装置118は、例えば、キーボード、マウス、タッチパネル、音声認識装置などであるが、これらに限られない。出力回路117は、管理サーバ101からのデータ出力を実現する。出力回路117は、外部ポートとしてHDMI(登録商標)またはDisplayPort等を含む。例えば、表示装置119は、出力回路117に接続される。表示装置119は、例えば、LCD(液晶ディスプレイ)、有機EL(有機エレクトロルミネッセンス)ディスプレイまたはプロジェクタ等であるが、これらに限られない。
管理者は、操作装置118および表示装置119を使って、管理サーバ101のメンテナンスを行うことができる。なお、操作装置118および表示装置119は、管理サーバ101に組み込まれたものであってもよい。また、操作装置118および表示装置119は、管理サーバ101に接続されていなくてもよい。例えば、管理者は、ネットワーク102と通信可能な端末装置を用いて管理サーバ101のメンテナンスを行ってもよい。
図10は、管理サーバ101の記憶部114に記憶されるデータを示す図である。記憶部114は、例えば、問題データ114Aと、計算データ114Bと、管理プログラム114Cと、変換プログラム114Dと、制御プログラム114Eとを記憶する。例えば、問題データ114Aは、組合せ最適化問題のデータを含む。例えば、計算データ114Bは、各計算サーバ103から収集された計算結果を含む。例えば、管理プログラム114Cは、管理部111の機能を実現するプログラムである。例えば、変換プログラム114Dは、変換部112の機能を実現するプログラムである。例えば、制御プログラム114Eは、制御部113の機能を実現するプログラムである。
図11は、計算サーバ103aの構成を示す図である。他の計算サーバ103は、計算サーバ103aと同様の構成であってもよいし、計算サーバ103aと異なる構成であってもよい。
計算サーバ103aは、例えば、通信回路131と、共有メモリ132と、プロセッサ133a~133dと、ストレージ134と、ホストバスアダプタ135とを備える。通信回路131、共有メモリ132、プロセッサ133a~133d、ストレージ134、ホストバスアダプタ135は、バス136を介して互いに接続される。
通信回路131は、ネットワーク102に接続された各装置との間でデータの送受信を行う。通信回路131は、例えば、有線LANのNIC(Network Interface Card)である。ただし、通信回路131は、無線LANなど、その他の種類の通信回路であってもよい。共有メモリ132は、プロセッサ133a~133dからアクセス可能なメモリである。共有メモリ132は、例えば、DRAMおよびSRAM等の揮発性メモリである。共有メモリ132は、不揮発性メモリ等の他の種類のメモリを含んでもよい。プロセッサ133a~133dは、共有メモリ132を介してデータを共有する。なお、共有メモリ132は、計算サーバ103aの全てメモリにより構成されていなくてもよい。例えば、計算サーバ103aの一部のメモリは、プロセッサ133a~133dのうちのいずれかからのみからアクセスできるローカルメモリであってもよい。
プロセッサ133a~133dは、計算処理を実行する電子回路である。プロセッサ133a~133dは、例えば、CPU、GPU、FPGA、ASICのいずれであってもよいし、これらの組合せであってもよい。また、プロセッサ133a~133dは、CPUコアまたはCPUスレッドであってもよい。プロセッサ133a~133dがCPUである場合、計算サーバ103aが備えるソケット数については、特に問わない。また、プロセッサ133a~133dは、PCI expressなどのバスを介して計算サーバ103aのその他の構成要素に接続されていてもよい。
図11の例では、計算サーバ103aは、4つのプロセッサ133a~133dを備える。しかし、1台の計算サーバ103aに含まれるプロセッサ数は、4個に限られない。
ストレージ134は、計算サーバ103aのプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、OSとアプリケーションの両方を含むものとする。ストレージ134は、揮発性メモリおよび不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリは、たとえば、DRAMまたはSRAMである。不揮発性メモリは、例えば、NANDフラッシュメモリ、NORフラッシュメモリ、ReRAMまたはMRAM等である。また、ストレージ134は、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置を含んでもよい。
ホストバスアダプタ135は、他の計算サーバ103との間のデータ通信を実現する。ホストバスアダプタ135は、ケーブル104aを介してスイッチ105に接続されている。ホストバスアダプタ135は、例えば、HCA(Host Channel Adaptor)である。ホストバスアダプタ135、ケーブル104a、スイッチ105で高スループットを実現可能なインターコネクトを形成することにより、並列的な計算処理の速度を向上させることができる。
図12は、ストレージ134に記憶されるデータを示す図である。ストレージ134は、例えば、計算データ134Aと、計算プログラム134Bと、制御プログラム134Cとを記憶する。計算データ134Aは、計算サーバ103aの計算途中のデータまたは計算結果を含む。計算データ134Aの少なくとも一部は、共有メモリ132、プロセッサのキャッシュまたはプロセッサのレジスタ等の、異なる記憶階層に記憶されてもよい。計算プログラム134Bは、各プロセッサ133における計算処理および、共有メモリ132およびストレージ134へのデータの保存処理を実現するプログラムである。制御プログラム134Cは、管理サーバ101の制御部113から送信された指令に基づき、計算サーバ103aを制御し、計算サーバ103aの計算結果を管理サーバ101に送信するプログラムである。
このような計算サーバ103aは、プロセッサ133a~133dが組合せ最適化問題を解くためのプログラムを実行する。このプログラムは、計算サーバ103aを、入力部12、更新部14および出力部16として機能させる。
このような構成の情報処理システム100は、PCクラスタとして利用することが可能である。PCクラスタは、複数台のコンピュータを接続し、1台のコンピュータでは得られない計算性能を実現するシステムである。情報処理システム100は、例えばMPI(Message Passing Interface)を使うことにより、複数の計算サーバ103にメモリが分散して配置されている構成であっても、並列的な計算を実行することが可能である。
なお、情報処理システム100は、高速リンクで接続された複数のGPUであってもよい。この場合、複数のGPUのそれぞれは、計算サーバ103と同様の処理を実行する。
本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組合せにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組合せてもよい。
10 求解装置
12 入力部
14 更新部
16 出力部
12 入力部
14 更新部
16 出力部
Claims (16)
- 0-1組合せ最適化問題を、前記0-1組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される1または複数の制約条件の下で解く求解装置であって、
第1変数および第2変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第1変数および前記第2変数を交互に更新する更新部と、
前記終了時刻における前記複数の要素のそれぞれの前記第1変数に基づき前記0-1組合せ最適化問題の解を出力する出力部と、
を備え、
前記複数の要素は、前記複数の離散変数に対応し、
前記第1変数および前記第2変数のそれぞれは、実数により表され、
前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、前記第1変数を前記第2変数に基づき更新し、
前記複数の要素のそれぞれについて、前記第2変数を前記第1変数に基づき更新し、
前記1または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第1変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第2変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
求解装置。 - 前記出力部は、
前記終了時刻における前記複数の要素のそれぞれについて、前記第1変数を予め設定されたしきい値により2値化した離散変数の値を算出し、
算出した前記複数の離散変数の値を前記0-1組合せ最適化問題の解として出力する
請求項1に記載の求解装置。 - 前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第1変数を更新した後に前記第2変数を更新し、
前記第2変数を更新した後に前記第2変数から前記補正値を減算する
請求項2に記載の求解装置。 - 前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第2変数を更新した後に前記第1変数を更新し、
前記第2変数を更新した後、且つ、前記第1変数を更新する前に、前記第2変数から前記補正値を減算する
請求項2に記載の求解装置。 - 前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、対象時刻における前記第1変数を、前記対象時刻より単位時間前の直前時刻における前記第1変数に、前記第2変数と予め定められた定数と前記単位時間とを乗算した値を加算することにより、算出する
請求項2に記載の求解装置。 - 前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、
前記複数の要素のそれぞれの前記第1変数と、対象要素と前記複数の要素のそれぞれとの組毎に前記0-1組合せ最適化問題により定められる作用係数とに基づき、外力を算出し、
時間経過に従って増加する関数に基づき定まる値と、前記対象要素の前記第1変数とを乗算した時間発展値を算出し、
前記直前時刻における前記第2変数に、前記時間発展値と前記外力とを加算した値に前記単位時間を乗算した値を、加算することにより、前記対象時刻における前記第2変数を算出する
請求項5に記載の求解装置。 - 前記0-1組合せ最適化問題は、N個の離散変数を含み、
前記更新部は、前記N個の離散変数のうちのi番目の離散変数に対応するi番目の要素の、前記対象時刻における前記第1変数を式(101)または式(102)により算出し、
iは、1からNまでの整数であり、
Dは、前記定数であり、
Δtは、前記単位時間であり、
tは、前記直前時刻であり、
t+Δtは、前記対象時刻であり、
xi(t)は、前記i番目の要素の前記直前時刻における前記第1変数であり、
yi(t)は、前記i番目の要素の前記直前時刻における前記第2変数であり、
xi(t+Δt)は、前記i番目の要素の前記対象時刻における前記第1変数であり、
yi(t+Δt)は、前記i番目の要素の前記対象時刻における前記第2変数である
請求項6に記載の求解装置。 - 前記0-1組合せ最適化問題は、0-1二次計画問題であり、
前記更新部は、前記i番目の要素の前記対象時刻における前記第2変数を式(103)または式(104)により算出し、
fi(t+Δt)は、式(105)により表され、fi(t)は、式(106)により表され、
zi(t)は、式(108)により表され、
hiは、N個の局所磁場係数を含む予め定められた配列に含まれる、i番目の局所磁場係数であり、
Ji,jは、N×N個の結合係数を含む予め定められた行列に含まれる、i行、j列の結合係数であり、
cは、予め定められた実数であり、
xj(t)は、前記N個の離散変数のうちのj番目の離散変数に対応するj番目の要素の前記直前時刻における前記第1変数であり、
xj(t+Δt)は、前記j番目の要素の前記対象時刻における前記第1変数であり、
α(t)は、tを変数とする予め定められた関数である
請求項7に記載の求解装置。 - 情報処理装置により、0-1組合せ最適化問題を、前記0-1組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される1または複数の制約条件の下で解く求解方法であって、
前記情報処理装置により、第1変数および第2変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第1変数および前記第2変数を交互に更新する更新ステップと、
前記情報処理装置により、前記終了時刻における前記複数の要素のそれぞれの前記第1変数に基づき前記0-1組合せ最適化問題の解を出力する出力ステップと、
を含み、
前記複数の要素は、前記複数の離散変数に対応し、
前記第1変数および前記第2変数のそれぞれは、実数により表され、
前記単位時間毎の更新ステップにおいて、前記情報処理装置は、
前記複数の要素のそれぞれについて、前記第1変数を前記第2変数に基づき更新し、
前記複数の要素のそれぞれについて、前記第2変数を前記第1変数に基づき更新し、
前記1または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第1変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第2変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
求解方法。 - 情報処理装置を、0-1組合せ最適化問題を、前記0-1組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される1または複数の制約条件の下で解く求解装置として機能させるためのプログラムであって、
前記情報処理装置を、
第1変数および第2変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第1変数および前記第2変数を交互に更新する更新部と、
前記終了時刻における前記複数の要素のそれぞれの前記第1変数に基づき前記0-1組合せ最適化問題の解を出力する出力部と、
して機能させ、
前記複数の要素は、前記複数の離散変数に対応し、
前記第1変数および前記第2変数のそれぞれは、実数により表され、
前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、前記第1変数を前記第2変数に基づき更新し、
前記複数の要素のそれぞれについて、前記第2変数を前記第1変数に基づき更新し、
前記1または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第1変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第2変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
プログラム。
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CA3219868A CA3219868A1 (en) | 2021-05-24 | 2022-04-14 | Solution finding device, solution finding method, and program |
US18/517,147 US20240095300A1 (en) | 2021-05-24 | 2023-11-22 | Solution finding device, solution finding method, and computer program product |
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2021087125A JP2022180174A (ja) | 2021-05-24 | 2021-05-24 | 求解装置、求解方法およびプログラム |
JP2021-087125 | 2021-05-24 |
Related Child Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
US18/517,147 Continuation US20240095300A1 (en) | 2021-05-24 | 2023-11-22 | Solution finding device, solution finding method, and computer program product |
Publications (1)
Publication Number | Publication Date |
---|---|
WO2022249785A1 true WO2022249785A1 (ja) | 2022-12-01 |
Family
ID=84229913
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
PCT/JP2022/017778 WO2022249785A1 (ja) | 2021-05-24 | 2022-04-14 | 求解装置、求解方法およびプログラム |
Country Status (4)
Country | Link |
---|---|
US (1) | US20240095300A1 (ja) |
JP (1) | JP2022180174A (ja) |
CA (1) | CA3219868A1 (ja) |
WO (1) | WO2022249785A1 (ja) |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2020196862A1 (ja) * | 2019-03-28 | 2020-10-01 | 株式会社 東芝 | 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム |
-
2021
- 2021-05-24 JP JP2021087125A patent/JP2022180174A/ja active Pending
-
2022
- 2022-04-14 CA CA3219868A patent/CA3219868A1/en active Pending
- 2022-04-14 WO PCT/JP2022/017778 patent/WO2022249785A1/ja active Application Filing
-
2023
- 2023-11-22 US US18/517,147 patent/US20240095300A1/en active Pending
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2020196862A1 (ja) * | 2019-03-28 | 2020-10-01 | 株式会社 東芝 | 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム |
Non-Patent Citations (1)
Title |
---|
GOTO, HAYATO: "Simulated branch algorithm", ABSTRACTS OF THE 2017 ANNUAL SPRING MEETING OF THE OPERATIONS RESEARCH SOCIETY OF JAPAN; MARCH 15TH - MARCH 17TH, 2017, THE OPERATIONS RESEARCH SOCIETY OF JAPAN, JP, 1 January 2020 (2020-01-01) - 13 March 2020 (2020-03-13), JP, pages 298 - 299, XP009541273, ISSN: 1883-1893 * |
Also Published As
Publication number | Publication date |
---|---|
US20240095300A1 (en) | 2024-03-21 |
CA3219868A1 (en) | 2022-12-01 |
JP2022180174A (ja) | 2022-12-06 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US20210224447A1 (en) | Grouping of pauli strings using entangled measurements | |
WO2020196862A1 (ja) | 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム | |
CN113379058B (zh) | 量子模拟方法及装置、电子设备及存储介质 | |
JP7297540B2 (ja) | 情報処理装置、puboソルバ、情報処理方法およびプログラム | |
KR20210134724A (ko) | 시뮬레이션된 양자 알고리즘에 기초한 데이터 검색 방법 및 장치 및 디바이스 | |
WO2021220445A1 (ja) | 演算システム、情報処理装置、および最適解探索処理方法 | |
US20220012307A1 (en) | Information processing device, information processing system, information processing method, and storage medium | |
WO2022249785A1 (ja) | 求解装置、求解方法およびプログラム | |
US11966450B2 (en) | Calculation device, calculation method, and computer program product | |
US20220012306A1 (en) | Information processing device, information processing system, information processing method, and storage medium | |
JP7472062B2 (ja) | 計算装置、計算方法およびプログラム | |
JP7474242B2 (ja) | 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム | |
JP7421545B2 (ja) | 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム | |
CN117669751A (zh) | 量子电路模拟方法、装置及电子设备 | |
CN117313877A (zh) | 量子电路处理方法、装置及电子设备 | |
JP2023166001A (ja) | 計算装置、計算プログラム、記録媒体及び計算方法 | |
CN117313879A (zh) | 量子电路处理方法、装置及电子设备 | |
CN117521829A (zh) | 量子电路模拟方法、装置及电子设备 | |
CN117313881A (zh) | 量子电路的分类方法、装置及电子设备 | |
CN117313878A (zh) | 量子电路处理方法、装置及电子设备 | |
CN117313882A (zh) | 量子电路处理方法、装置及电子设备 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 22811059 Country of ref document: EP Kind code of ref document: A1 |
|
WWE | Wipo information: entry into national phase |
Ref document number: 3219868 Country of ref document: CA |
|
NENP | Non-entry into the national phase |
Ref country code: DE |
|
122 | Ep: pct application non-entry in european phase |
Ref document number: 22811059 Country of ref document: EP Kind code of ref document: A1 |