WO2022249785A1 - 求解装置、求解方法およびプログラム - Google Patents

求解装置、求解方法およびプログラム Download PDF

Info

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
Application number
PCT/JP2022/017778
Other languages
English (en)
French (fr)
Inventor
賢 鈴木
Original Assignee
株式会社東芝
東芝デジタルソリューションズ株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社東芝, 東芝デジタルソリューションズ株式会社 filed Critical 株式会社東芝
Priority to CA3219868A priority Critical patent/CA3219868A1/en
Publication of WO2022249785A1 publication Critical patent/WO2022249785A1/ja
Priority to US18/517,147 priority patent/US20240095300A1/en

Links

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)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational 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

不等式の制約条件の下で０－１組合せ最適化問題を解く。求解装置は、更新部と、出力部とを備える。更新部は、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、第１変数および第２変数を交互に更新する。出力部は、終了時刻における複数の要素のそれぞれの第１変数に基づき０－１組合せ最適化問題の解を出力する。単位時間毎の更新処理において、更新部は、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数を代入した不等式を満たさない場合、複数の要素のそれぞれの第２変数から、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。

Description

求解装置、求解方法およびプログラム
　本発明の実施形態は、求解装置、求解方法およびプログラムに関する。
　従来、一次不等式による制約の下でイジング問題を解く方法として、スラック変数を利用した定式化をする方法が知られている。この方法は、一次不等式の定数項をＷｍとした場合、目的関数に、Ｗ個（Ｗは自然数）のスラック変数により表されたペナルティ項を加算して、イジング問題を解く。従って、この方法は、目的関数にＮ個の決定変数が含まれるならば、Ｎ＋Ｗ個のイジングスピンを含むイジング問題を解く。また、例えば、Ｍ個の一次不等式による制限の下でイジング問題を解くのであれば、この方法は、Ｎ＋｛Ｗ１＋…＋Ｍｍ＋…＋ＷＭ｝個のイジングスピンを含むイジング問題を解く。なお、Ｍｍは、Ｍ個の一次不等式のうちのｍ番目の一次不等式の定数項であり、自然数を表す。
　株式のポートフォリオ最適化問題は、リスク回避のために、銘柄の投資割合を決定変数とし、複数の決定変数のそれぞれに重みを乗算した合計を上限値以下とする制約条件の下で、解かれる場合が多い。制約条件は、下記の式のような一次不等式により表される。
　ｘＡ１・ωＡ１＋ｘＡ２・ωＡ２＋…＋ｘＳｍ・ωＳｍ＋…≦Ｗｍ
　この式において、Ｗｍは、上限値を表す。ｘＡ１は、銘柄Ａの投資割合を例えば１％とする決定変数を表す。ｘＡ２は、銘柄Ａの投資割合を例えば２％とする決定変数を表す。ｘＳｓは、銘柄Ｓの投資割合をｓ％とする決定変数を表す。また、ωＡ１、ωＡ２およびωＳｍは、重みを表す。
　例えば、株式のポートフォリオ最適化問題は、複数の銘柄を、業種または企業規模毎にグループ化し、グループ毎に上記の制約条件を設定する。例えば、東京証券取引所は、約２０００銘柄を取り扱っている。また、東京証券取引所は、これら約２０００銘柄を、３３種類の業種および３種類の規模に分類している。従って、東京証券所で取り扱われている約２０００銘柄は、３３業種×３規模の合計９９個のグループに分けることができる。
　１つの銘柄の投資割合を１％から１００％までの１００段階で離散化して、東京証券取引所により取り扱われる２０００銘柄に対してポートフォリオ最適化問題を解く場合、目的関数に含まれる決定変数の個数は、２０００００個となる。
　これに対して、９９個のグループの全てに対して一次不等式による制約条件を設定するとする。上限値（Ｗｍ）を１０００段階で離散化して、スラッグ変数を利用した定式化をした場合、スラッグ変数の個数は、９９０００個となる。このため、９９個のグループの全てに対して一次不等式による制約条件を設定したポートフォリオ最適化問題を解く場合、イジングスピンの個数が２００００＋９９０００＝２９９０００個のイジング問題を解かなければならない。
　このように、スラック変数を利用した定式化をしてポートフォリオ最適化問題を解く場合、制約条件が無い場合と比較して、例えば東京証券取引所の例では、イジングスピンの数が５０％程度増加してしまう。ソルバー等によりイジング問題を解く場合、イジングスピンの数が増加してしまうと、演算量および演算時間が多くなり、求解のためのコストが大きくなってしまう。
特開２０１７－７３１０６号公報 特開２０１９－１４５０１０号公報 特開２０２１－０４３５８９号公報 国際公開第２０２０／１９６９１５号
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］，［２０２１年５月１７日検索］，インターネット，　<http://www.ynl.t.u-tokyo.ac.jp/project/RobotBrainCREST/publications/pdf/tsuchiya/4_01.pdf> 土屋和雄，西山岳宏，辻田勝吉，　"分岐特性を用いた組合せ最適化問題の近似解法　第２報：決定論的アニーリングアルゴリズムの解析"，　［online］，［２０２１年５月１７日検索］，インターネット，　<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
　本発明が解決しようとする課題は、不等式により表される１または複数の制約条件の下で０－１組合せ最適化問題を少ない演算コストで解くことができる求解装置、求解方法およびプログラムを提供する。
　実施形態に係る求解装置は、０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く。前記求解装置は、更新部と、出力部とを備える。前記更新部は、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する。前記出力部は、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する。前記複数の要素は、前記複数の離散変数に対応する。前記第１変数および前記第２変数のそれぞれは、実数により表される。前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新する。前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新する。前記単位時間毎の更新処理において、前記更新部は、前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。
図１は、シミュレーテッド分岐アルゴリズムの分岐現象を示す図である。 図２は、本実施形態に係る求解装置の機能構成を示す図である。 図３は、更新部の処理の流れの第１例を示すフローチャートである。 図４は、更新部の処理の流れの第２例を示すフローチャートである。 図５は、制約処理の第１例の処理の流れを示すフローチャートである。 図６は、制約処理の第２例の処理の流れを示すフローチャートである。 図７は、制約処理の第３例の処理の流れを示すフローチャートである。 図８は、情報処理システムの構成図である。 図９は、管理サーバの構成図である。 図１０は、記憶部に記憶されるデータを示す図である。 図１１は、計算サーバの構成図である。 図１２は、ストレージに記憶されるデータを示す図である。
　（０－１組合せ最適化問題）
　イジング問題を解くために使われる装置の一例として、イジングマシンが挙げられる。イジングマシンは、イジングモデルの基底状態のエネルギーを計算する。これまで、イジングモデルは、主に強磁性体や相転移現象のモデルとして使われることが多かった。しかし、近年、イジングモデルは、０－１組合せ最適化問題を解くためのモデルとしての利用が増えている。式（１）は、イジングモデルのエネルギーを示す。
Figure JPOXMLDOC01-appb-M000015
　ｓｉ、ｓｊはスピンを表す。スピンは、＋１または－１のいずれかの値をとる２値変数である。ｓｉは、ｉ番目のスピンを表す。ｓｊは、ｊ番目のスピンを表す。ｉおよびｊは、１以上、Ｎ以下の整数である。Ｎは、スピンの数を表し、２以上の整数である。ｈｉは、ｉ番目のスピンに作用する局所磁場を表す。Ｊは、２つのスピン間に作用する力を表す結合係数の行列である。Ｊは、対角成分が０である実対称行列である。Ｊｉｊは、Ｊのｉ行ｊ列の要素を表す。つまり、Ｊｉｊは、ｉ番目のスピンと、ｊ番目のスピンとの間に作用する力を表す結合係数である。
　イジングマシンは、式（１）により表されるエネルギーＥＩｓｉｎｇを目的関数とし、エネルギーＥＩｓｉｎｇを可能な限り小さくする解を算出する。エネルギーＥＩｓｉｎｇが最小値となるイジングモデルの解（ｓ１、ｓ２、・・・、ｓＮ）は、最適解と呼ばれる。ただし、イジングモデルの解は、最適解ではなく、エネルギーＥＩｓｉｎｇが最小値に近い近似解であってもよい。すなわち、イジング問題は、最適解のみならず、近似解を算出する問題であってもよい。
　また、０または１のいずれかの値をとる離散変数（ビット）の２次関数を目的関数とする０－１組合せ最適化問題は、０－１二次計画問題と呼ばれる。離散変数（ビット）は、（１＋ｓｉ）／２の演算を用いることにより、ｓｉに変換される。つまり、０－１二次計画問題は、式（１）で表されるイジング問題と等価であるといえる。従って、０－１二次計画問題は、イジング問題に変換し、イジングマシンにより解を算出することが可能である。
　イジングマシンは、例えば、量子アニーラ、コヒーレントイジングマシンおよび量子分岐マシン等によりハードウェア実装される。量子アニーラは、超伝導回路を使って量子アニーリングを実現する。コヒーレントイジングマシンは、光パラメトリック発振器で形成されたネットワークの発振現象を利用する。量子分岐マシンは、カー効果を有するパラメトリック発振器のネットワークにおける量子力学的な分岐現象を利用する。これらのハードウェア実装されたイジングマシンは、計算時間の大幅な短縮を実現する可能性がある一方、大規模化および安定的な運用が難しいという課題もある。
　イジング問題は、広く普及しているデジタルコンピュータを用いて解を算出することも可能である。デジタルコンピュータは、量子アニーラ、コヒーレントイジングマシンおよび量子分岐マシン等と比べ、大規模化および安定運用が可能である。シミュレーテッドアニーリング（ＳＡ）は、デジタルコンピュータでイジング問題を解くためのアルゴリズムの一例である。ただし、シミュレーテッドアニーリングは、それぞれの変数が逐次更新される逐次更新アルゴリズムであるため、並列化による計算処理の高速化は難しい。
　非特許文献２には、０－１組合せ最適化問題を解くためのアルゴリズムとして、シミュレーテッド分岐アルゴリズムが提案されている。シミュレーテッド分岐アルゴリズムは、イジングモデルを用いて、デジタルコンピュータによって、規模の大きい０－１組合せ最適化問題を高速に解くことが可能である。シミュレーテッド分岐アルゴリズムは、ＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）、マイクロプロセッサ、ＧＰＵ（Ｇｒａｐｈｉｃｓ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）、または、これらの組合せの回路等の電子回路によっても、規模の大きい０－１組合せ最適化問題を高速に解くことが可能である。
　（シミュレーテッド分岐アルゴリズム）
　シミュレーテッド分岐アルゴリズムは、それぞれがＮ個の要素に対応する変数ｘｉおよび変数ｙｉを用いる。変数ｘｉを第１変数、変数ｙｉを第２変数と呼ぶ場合もある。シミュレーテッド分岐アルゴリズムにおいて、Ｎ個の要素のそれぞれは、仮想的な粒子を表す。Ｎ個の要素は、イジング問題のＮ個のスピンに対応する。従って、Ｎ個の要素は、組合せ最適化問題のＮ個の離散変数（ビット）に対応する。変数ｘｉおよび変数ｙｉは、いずれも、実数で表される連続変数である。変数ｘｉは、Ｎ個の粒子のうちのｉ番目の粒子の位置を表す。変数ｙｉは、ｉ番目の粒子の運動量を表す。Ｎは、２以上の整数である。ｉは、１以上、Ｎ以下の整数を表し、Ｎ個の要素のそれぞれを特定するインデックスを表す。
　シミュレーテッド分岐アルゴリズムは、それぞれＮ個ある変数ｘｉおよび変数ｙｉについて、下記の式（２）の連立常微分方程式を数値的に解く。
Figure JPOXMLDOC01-appb-M000016
　Ｈは、下記の式（３）のハミルトニアンである。
Figure JPOXMLDOC01-appb-M000017
　係数Ｄは、予め定められた定数であり、離調（ｄｅｔｕｎｉｎｇ）に相当する。係数ｐ（ｔ）は、ポンピング振幅（ｐｕｍｐｉｎｇ　ａｍｐｌｉｔｕｄｅ）に相当し、シミュレーテッド分岐アルゴリズムの計算時に更新回数に応じて値が単調増加する。ｔは、時刻を表す変数である。係数ｐ（ｔ）の初期値は０に設定されていてもよい。係数Ｋは、予め定められた定数であって、正のカー係数（Ｋｅｒｒ　ｃｏｅｆｆｉｃｉｅｎｔ）に相当する。なお、係数Ｋは、０であってもよい。
　ｆｉは、外力を表し、下記の式（４）で表される。
Figure JPOXMLDOC01-appb-M000018
　式（４）のｚｉは、式（３）の中の小カッコの内の数式を変数ｘｉで偏微分した式である。式（３）の中の小カッコの内の数式は、イジングモデルのエネルギーＥＩｓｉｎｇに対応する。
　ｃは、係数である。ｃは、例えば、計算を実行する前に予め定められる定数であってもよい。また、α（ｔ）は、ｐ（ｔ）とともに増加する係数である。
　そして、シミュレーテッド分岐アルゴリズムは、ｐ（ｔ）の値を初期値（例えば、０）から所定の値まで増加させた後における変数ｘｉの符号に基づき、スピンｓｉの値を算出する。シミュレーテッド分岐アルゴリズムは、例えば、ｘｉ＞０の場合にｓｇｎ（ｘｉ）＝１、ｘｉ＜０の場合にｓｇｎ（ｘｉ）＝－１となる符号関数を用いて、スピンｓｉの値を算出する。
　（シミュレーテッド分岐アルゴリズムの演算）
　シミュレーテッド分岐アルゴリズムは、シンプレクティック・オイラー法を用いて、式（２）、式（３）および式（４）によって与えられる微分方程式を解く。
　ここで、シンプレクティック・オイラー法を使う場合、式（２）、式（３）および式（４）によって与えられる微分方程式は、式（５）または式（６）に示すような、離散的な漸化式に書き換えられる。
Figure JPOXMLDOC01-appb-M000019
Figure JPOXMLDOC01-appb-M000020
　ｔは、時刻を表す。Δｔは、単位時間（時間ステップ、時間刻み幅）を表す。
　シミュレーテッド分岐アルゴリズムを実行する場合、デジタルコンピュータまたはＦＰＧＡ等の電子回路は、式（５）または式（６）のアルゴリズムに基づき、それぞれＮ個ある変数ｘｉおよび変数ｙｉを初期時刻から単位時間毎に順次に、且つ、変数ｘｉと変数ｙｉとを交互に、更新する。そして、デジタルコンピュータまたはＦＰＧＡ等の電子回路は、終了時刻におけるＮ個の変数ｘｉの値を、符号関数を用いて２値化して、Ｎ個のスピンの値を出力する。
　なお、式（５）および式（６）は、微分方程式との対応関係を示すために、時刻ｔおよび単位時間Δｔを用いて表されている。ただし、シンプレクティック・オイラー法をデジタルコンピュータまたはＦＰＧＡ等の電子回路で実行する場合、式（５）および式（６）を演算するためのアルゴリズムは、明示的なパラメータとして時刻ｔおよび単位時間Δｔを含まなくてよい。例えば、単位時間Δｔを１とする場合、式（５）および式（６）を演算するためのアルゴリズムは、単位時間Δｔを含まなくてよい。例えば、明示的なパラメータとして時刻ｔを含まない場合、式（５）および式（６）を演算するアルゴリズムは、ｘｉ（ｔ＋Δｔ）をｘｉ（ｔ）の更新後の値として処理を実行する。すなわち、式（５）および式（６）を演算するアルゴリズムは、“ｔ”を更新前の変数を特定するパラメータ、“ｔ＋Δｔ”を更新後の変数を特定するパラメータとして処理を実行する。以降で説明する式（５）および式（６）を改良したアルゴリズムも同様である。
　図１は、シミュレーテッド分岐アルゴリズムにより最適化問題を解いた場合における、変数ｘｉの分岐現象を表す図である。シミュレーテッド分岐アルゴリズムにより最適化問題を解いた場合、系のパラメータが変化することに伴い、安定運動状態が１個のみの系から、安定運動状態が２個の系へと遷移する分岐現象が生じる。図１に示すように、分岐現象が進むと、変数ｘｉは、－１または＋１の近傍に集中するが、－１より小さい領域、または、＋１より大きい領域にも広がる。
　（改良アルゴリズム）
　発明者は、上述のシミュレーテッド分岐アルゴリズムを改良して、不等式により表される１または複数の制約条件の下で、０－１組合せ最適化問題を解く改良アルゴリズムを発明した。
　改良アルゴリズムは、式（３）のハミルトニアンに代えて、式（７）に示すハミルトニアンを用いて、式（２）の連立常微分方程式を数値的に解く。
Figure JPOXMLDOC01-appb-M000021
　Ｇｍは、１または複数の制約条件のうちのｍ番目の制約条件に対応するペナルティ項を表す。１または複数の制約条件のそれぞれは、目的関数に含まれる複数の離散変数を含む不等式により表される。例えば、不等式は、複数の離散変数の一次項を含む一次不等式、または、複数の離散変数の二次項を含む二次不等式である。
　ペナルティ項は、対応する不等式を満たす場合に、ハミルトニアンにエネルギーとして０を加え、対応する不等式を満たさない場合に、ハミルトニアンに正のエネルギーを加える。より具体的には、Ｇｍは、制約条件を表す不等式の境界から、Ｎ個の要素により特定される位置までの距離が大きい程、大きなエネルギーをハミルトニアンに加える関数である。なお、不等式の境界は、不等式を満たすか満たさないかの境界である。すなわち、不等式の境界は、不等式における不等号を、等号により代えた方程式により表される。
　そして、改良アルゴリズムは、シンプレクティック・オイラー法により、連立常微分方程式を式（８）または式（９）に示すような離散的な漸化式を用いて解く。
Figure JPOXMLDOC01-appb-M000022
Figure JPOXMLDOC01-appb-M000023
　式（８）および式（９）に示す漸化式は、式（５）および式（６）に示す漸化式と比較すると、第２変数（ｙｉ（ｔ＋Δｔ））の算出式が異なる。具体的には、式（８）および式（９）に示す漸化式は、第２変数（ｙｉ（ｔ＋Δｔ））の算出式に、１または複数の不等式のそれぞれのペナルティ項（Ｇｍ）のｘｉについての偏微分値の合計に、Δｔを乗算した値が減算されている点において異なる。
　１または複数の不等式のそれぞれの偏微分値は、対応する不等式の境界からＮ個の要素により特定される位置までの距離における対応する要素の成分に相当する。この偏微分値は、対応する不等式の境界からＮ個の要素により特定される位置までの距離における対応する要素の成分が大きい程、大きい値となる。
　ただし、ペナルティ項は、対応する不等式を満たす場合に、０となる。従って、１または複数の不等式のそれぞれの偏微分値は、Ｎ個の要素が不等式を満たしている場合には、０となる。
　ここで、シミュレーテッド分岐アルゴリズムは、時刻を表す変数であるｔを、単位時間（Δｔ）毎に増加させながら、複数の第１変数（ｘｉ）および複数の第２変数（ｙｉ）を交互に更新する。従って、シミュレーテッド分岐アルゴリズムは、単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれについて、不等式を満たすか否かを判断することができる。
　そこで、改良アルゴリズムは、単位時間（Δｔ）毎に、複数の第２変数（ｙｉ）を更新した後に、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数（ｘｉ）を代入した不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数を代入した不等式を満たさない場合、複数の第２変数（ｙｉ）のそれぞれから、対応する補正値を減算する。つまり、改良アルゴリズムは、単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数を代入した不等式を満たさない場合、複数の第２変数（ｙｉ）のそれぞれから、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。
　これにより、改良アルゴリズムは、複数の第１変数（ｘｉ）および複数の第２変数（ｙｉ）を単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれを満たすように、複数の第１変数（ｘｉ）および複数の第２変数（ｙｉ）を補正することができる。この結果、改良アルゴリズムは、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。
　このような改良アルゴリズムは、不等式により表される１または複数の制約条件の下で、０－１組合せ最適化問題を解く場合であっても、制約条件を設定しない場合と同一の変数の数により、解を得ることができる。従って、このような改良アルゴリズムは、スラッグ変数を設けずに、簡易に、不等式により表される１または複数の制約条件の下で０－１組合せ最適化問題を解くことができる。
　（制約条件が一次不等式であるのペナルティ項および偏微分値の第１例）
　１または複数の制約条件の一つは、一次不等式により表される。例えば、改良アルゴリズムは、Ｍ個（Ｍは１以上の整数）の一次不等式の制約の下で、０－１組合せ最適化問題を解く。この場合において、１または複数の制約条件のうちのｍ番目の制約条件は、式（１０）により表される。
Figure JPOXMLDOC01-appb-M000024
　式（１０）において、ｍは、１以上、Ｍ以下の整数である。式（１０）において、ｘｉは、ｉ番目の離散変数である。式（１０）において、ｘｊは、ｊ番目の離散変数である。ωｍ，ｉは、ｉ番目の離散変数に乗算される係数である。Ｗｍは、予め定められた定数である。
　制約条件が一次不等式により表される場合、ペナルティ項であるＧｍは、式（１１）に示すように表すことができる。
Figure JPOXMLDOC01-appb-M000025
　Ａｍは、ｍ番目の制約条件に対して予め定められた係数である。ｋは、１以上の整数である。例えば、ｋは、２であってもよい。
　Ｓｍは、一次不等式における定数項（Ｗｍ）を除く関数における複数の離散変数のそれぞれに、対応する第１変数（ｘｉ）を代入して得られる値である。具体的には、Ｓｍは、式（８）の漸化式を適用する場合には、式（１２）のように表される。また、Ｓｍは、式（９）の漸化式を適用する場合には、式（１３）のように表される。
Figure JPOXMLDOC01-appb-M000026
Figure JPOXMLDOC01-appb-M000027
　この場合、偏微分値は、式（１４）に示すように表される。
Figure JPOXMLDOC01-appb-M000028
　そこで、制約条件が一次不等式で表される場合、改良アルゴリズムは、単位時間（Δｔ）毎に、複数の第２変数（ｙｉ）を更新した後に、複数の離散変数のそれぞれに対応する第１変数（ｘｉ）を代入した一次不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間（Δｔ）毎に、対応する一次不等式を満たさない場合、複数の第２変数（ｙｉ）のそれぞれから、式（１４）に表される偏微分値に単位時間を表すΔｔを乗算した補正値を、減算する。
　このような処理を実行する改良アルゴリズムは、制約条件を設定しない場合と同一の変数の数により、一次不等式により表される制約条件の下で０－１組合せ最適化問題を簡易に解くことができる。
　（制約条件が一次不等式であるのペナルティ項および偏微分値の第２例）
　式（１４）に示す偏微分値は、Ａｍが大きい場合、計算誤差が蓄積して不安定となる可能性がある。そこで、改良アルゴリズムは、制約条件が一次不等式である場合、式（１５）により表される偏微分値を用いてもよい。すなわち、改良アルゴリズムは、単位時間（Δｔ）毎に、対応する一次不等式を満たさない場合、複数の第２変数（ｙｉ）のそれぞれから、式（１５）に表される偏微分値に単位時間を表すΔｔを乗算した補正値を、減算してもよい。
Figure JPOXMLDOC01-appb-M000029
　Ｂｍは、０以上１以下の値であり、ｍ番目の制約条件に対して予め定められた係数である。Ｅｍは、一次不等式における定数項（Ｗｍ）を除く関数における複数の離散変数のそれぞれに、対応する第２変数（ｙｉ）を代入して得られる値である。具体的には、Ｅｍは、式（８）または式（９）の漸化式を適用する場合には、式（１６）のように表される。
Figure JPOXMLDOC01-appb-M000030
　Ｅｍは、複数の要素により特定される方向の、一次不等式の境界に対する傾きを表す。従って、式（１５）の偏微分値を用いた場合、補正値は、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分と、複数の要素により特定される方向の、一次不等式の境界に対する傾きの対応する要素の成分とを加算した値に応じた値となる。これにより、改良アルゴリズムは、Ａｍが大きいために距離に応じた成分が不安定となる場合であっても、補正値を安定化させることができる。
　（制約条件が二次不等式であるのペナルティ項および偏微分値）
　１または複数の制約条件の一つは、二次不等式により表されてもよい。この場合において、制約条件は、式（１７）により表される。
Figure JPOXMLDOC01-appb-M000031
　式（１７）において、ｘｉは、ｉ番目の離散変数である。式（１７）において、ｘｊは、ｊ番目の離散変数である。Ｑｉ，ｊは、Ｎ行、Ｎ列の半正定値行列に含まれるｉ行、ｊ列の値である。ｑは、予め定められた定数である。
　制約条件が二次不等式により表される場合、ペナルティ項であるＧは、式（１８）に示すように表すことができる。
Figure JPOXMLDOC01-appb-M000032
　Ａは、二次不等式による制約条件に対して予め定められた係数である。ｋは、１以上の整数である。例えば、ｋは、２であってもよい。
　Ｓは、二次不等式における定数項（ｑ）を除く関数における複数の離散変数のそれぞれに、対応する第１変数（ｘｉ，ｘｊ）を代入して得られる値である。具体的には、Ｓは、式（８）の漸化式を適用する場合には、式（１９）のように表される。また、Ｓは、式（９）の漸化式を適用する場合には、式（２０）のように表される。
Figure JPOXMLDOC01-appb-M000033
Figure JPOXMLDOC01-appb-M000034
　この場合、偏微分値は、式（８）の漸化式を適用する場合には、式（２１）に示すように表される。また、偏微分値は、式（９）の漸化式を適用する場合には、式（２２）に示すように表される。
Figure JPOXMLDOC01-appb-M000035
Figure JPOXMLDOC01-appb-M000036
　そこで、制約条件が二次不等式で表される場合、改良アルゴリズムは、単位時間（Δｔ）毎に、複数の第２変数（ｙｉ）を更新した後に、複数の離散変数のそれぞれに対応する第１変数（ｘｉ，ｘｊ）を代入した二次不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間（Δｔ）毎に、対応する二次不等式を満たさない場合、複数の第２変数（ｙｉ）のそれぞれから、式（２１）または式（２２）に表される偏微分値に単位時間を表すΔｔを乗算した補正値を、減算する。
　イジングモデルは、二次不等式による制約条件の下で０－１組合せ最適化問題を解くことは困難である。これに対して、改良アルゴリズムは、二次不等式により表される制約条件の下で、０－１組合せ最適化問題を、簡易に解を得ることができる。
　なお、改良アルゴリズムは、１または複数の制約条件の中に、一次不等式による制約条件と、二次不等式による制約条件との両者が含まれていても、０－１組合せ最適化問題を解くことができる。また、改良アルゴリズムは、１または複数の制約条件に一次不等式による制約条件のみが含まれていてもよい。また、改良アルゴリズムは、１つの二次不等式による制約条件のみの下で、０－１組合せ最適化問題を解いてもよい。
　（機能ブロック構成）
　図２は、本実施形態に係る求解装置１０の機能構成を示す図である。
　求解装置１０は、改良アルゴリズムを用いて、０－１組合せ最適化問題を不等式により表される１または複数の制約条件の下で解く。求解装置１０は、コンピュータ等の情報処理装置、ネットワークを介して複数のコンピュータまたはサーバが相互に通信をして構成されるコンピュータシステム、または、複数台のコンピュータが連携して情報処理を実行するＰＣクラスタ等により実現される。また、求解装置１０は、ＣＰＵ、マイクロプロセッサ、ＧＰＵ、ＦＰＧＡまたはＡＳＩＣ、または、これらの組合せの回路等の電子回路によって実現される。
　求解装置１０は、機能構成として、入力部１２と、更新部１４と、出力部１６とを備える。
　入力部１２は、０－１組合せ最適化問題の目的関数を定義するための情報（例えば、Ｎ、Ｊ、ｈ）、および、改良アルゴリズムを実行するために必要な係数を表す情報（例えば、Ｄ、ｃ、Δｔ、Ｔ，ｐ（ｔ）、α（ｔ））を外部装置から受け取る。さらに、入力部１２は、１または複数の制約条件を定義するための情報（例えば、ωｍ，ｉ，Ｗｍ，Ｑｉ，ｊ，ｑ）、および、制約処理を実行するために必要な係数（例えば、Ａｍ，ｋ，Ｂｍ，Ａ）を表す情報を外部装置から受け取る。そして、入力部１２は、受け取ったこれらの情報を更新部１４に与える。
　更新部１４は、改良アルゴリズムを用いて、第１変数（ｘｉ）および第２変数（ｙｉ）が対応付けられた複数の要素のそれぞれについて、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎に順次に、第１変数（ｘｉ）および第２変数（ｙｉ）を交互に更新する。
　出力部１６は、終了時刻（ｔ＝Ｔ）における複数の要素のそれぞれの第１変数（ｘｉ）に基づき、０－１組合せ最適化問題の解を出力する。例えば、出力部１６は、終了時刻における複数の要素のそれぞれについて、第１変数（ｘｉ）を予め設定されたしきい値により２値化した離散変数の値を算出する。そして、出力部１６は、算出した複数の離散変数の値を０－１組合せ最適化問題の解として出力する。
　ここで、複数の要素は、０－１組合せ最適化問題の複数の離散変数に対応する。また、第１変数（ｘｉ）および第２変数（ｙｉ）のそれぞれは、実数により表される。
　そして、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第１変数（ｘｉ）を第２変数（ｙｉ）に基づき更新する。また、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第２変数（ｙｉ）を第１変数（ｘｉ）に基づき更新する。
　例えば、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第１変数（ｘｉ）を更新した後に第２変数（ｙｉ）を更新する。これに代えて、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第２変数（ｙｉ）を更新した後に第１変数（ｘｉ）を更新してもよい。
　さらに、単位時間毎の更新処理において、更新部１４は、第２変数（ｙｉ）を更新した後に、０－１組合せ最適化問題の解が１または複数の制約条件を満たすように、制約処理を実行する。更新部１４は、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数（ｘｉ）を代入した不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙｉ）から、対応する不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値（ｒｍ，ｉ）を減算する。
　単位時間毎の更新処理において、１または複数の制約条件のそれぞれについて制約処理を実行することにより、更新部１４は、複数の第１変数（ｘｉ）が制約条件を満たしていない場合、複数の第１変数（ｘｉ）が制約条件を満たす方向に変化するように、複数の要素に対応する複数の第２変数（ｙｉ）を補正することができる。これにより、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。
　なお、単位時間毎の更新処理において、更新部１４は、第１変数（ｘｉ）を更新した後に第２変数（ｙｉ）を更新する場合には、複数の要素のそれぞれについて、第２変数（ｙｉ）を更新した後に第２変数（ｙｉ）から補正値（ｒｍ，ｉ）を減算する。また、単位時間毎の更新処理において、更新部１４は、第２変数（ｙｉ）を更新した後に第１変数（ｘｉ）を更新する場合には、複数の要素のそれぞれについて、第２変数（ｙｉ）を更新した後、且つ、第１変数（ｘｉ）を更新する前に、第２変数（ｙｉ）から補正値（ｒｍ，ｉ）を減算する。
　（処理フロー）
　図３は、更新部１４の処理の流れの第１例を示すフローチャートである。更新部１４は、例えば、図３に示す流れで処理を実行する。
　まず、Ｓ１０１において、更新部１４は、０－１組合せ最適化問題を解くためのパラメータを設定する。具体的には、更新部１４は、Ｎ×Ｎ個の結合係数を含む行列であるＪ、および、Ｎ個の局所磁場を表す局所磁場係数を含む配列であるｈを設定する。さらに、更新部１４は、係数であるＤ、係数であるｃ、単位時間を表すΔｔ、終了時刻を表すＴ、関数であるｐ（ｔ）、および、関数であるα（ｔ）を設定する。ｐ（ｔ）およびα（ｔ）は、ｔ＝初期時刻（例えば０）で０、ｔ＝終了時刻（Ｔ）で１となる増加関数である。更新部１４は、Ｊ、ｈを入力部１２からの受け取った情報に応じて設定する。更新部１４は、Ｄ、ｃ、Δｔ、Ｔ、ｐ（ｔ）およびα（ｔ）を、入力部１２から受け取った値に応じて設定してもよいし、予め決定されており変更できない値を設定してもよい。
　続いて、Ｓ１０２において、更新部１４は、制約条件に関する情報を設定する。具体的には、更新部１４は、１または複数の制約条件を定義するための情報、および、制約処理を実行するために必要な係数を表す情報を設定する。例えば、更新部１４は、制約条件が一次不等式により表される場合、変数に乗算される係数であるωｍ，ｉ、定数項であるＷｍ、制約処理を実行するために必要な係数であるＡｍ（またはＡｍおよびＢｍ）、並びに、ｋを設定する。また、例えば、制約条件が二次不等式により表される場合、更新部１４は、変数に乗算される係数であるＱｉ，ｊおよび定数項であるｑ、制約処理を実行するために必要な係数であるＡおよびｋを設定する。なお、更新部１４は、制約処理を実行するために必要な係数であるＡｍ、Ｂｍ、ｋおよびＡを、入力部１２から受け取ったパラメータに応じて設定してもよいし、予め決定されており変更できない値を設定してもよい。
　続いて、Ｓ１０３において、更新部１４は、変数を初期化する。具体的には、更新部１４は、時刻を表す変数であるｔを初期時刻（例えば、０）に初期化する。さらに、更新部１４は、Ｎ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））のそれぞれおよびＮ個の第２変数（ｙ１（ｔ）～ｙＮ（ｔ））のそれぞれに、ユーザから受け取った初期値、予め定められた固定値、または、乱数を代入する。
　続いて、更新部１４は、Ｓ１０４とＳ１１６との間のループ処理を、ｔがＴより大きくなるまで繰り返す。１回のループ処理において、更新部１４は、対象時刻（ｔ＋Δｔ）におけるＮ個の第１変数（ｘ１（ｔ＋Δｔ）～ｘＮ（ｔ＋Δｔ））を、直前時刻（ｔ）におけるＮ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））、および、直前時刻（ｔ）におけるＮ個の第２変数（ｙ１（ｔ）～ｙＮ（ｔ））に基づき算出する。また、１回のループ処理において、更新部１４は、対象時刻（ｔ＋Δｔ）におけるＮ個の第２変数（ｙ１（ｔ＋Δｔ）～ｙＮ（ｔ＋Δｔ））を、対象時刻（ｔ＋Δｔ）におけるＮ個の第１変数（ｘ１（ｔ＋Δｔ）～ｘＮ（ｔ＋Δｔ））および直前時刻（ｔ）におけるＮ個の第２変数（ｙ１（ｔ）～ｙＮ（ｔ））に基づき算出する。
　なお、直前時刻（ｔ）は、対象時刻（ｔ＋Δｔ）より単位時間（Δｔ）前の時刻である。すなわち、更新部１４は、Ｓ１０４とＳ１１６との間のループ処理を繰り返すことにより、Ｎ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））およびＮ個の第２変数（ｙ１（ｔ）～ｙＮ（ｔ））を、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎に順次に更新する。
　続いて、更新部１４は、Ｓ１０５とＳ１０７との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。ｉは、１からＮまでの整数であり、Ｎ個の要素のうちの処理対象を表すインデックスである。Ｎ個の要素のそれぞれは、第１変数（ｘｉ（ｔ））および第２変数（ｙｉ（ｔ））が対応付けられる。Ｓ１０５とＳ１０７との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。
　Ｓ１０６において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘｉ（ｔ））に、対象要素の直前時刻（ｔ）における第２変数（ｙｉ（ｔ））と予め定められた定数（Ｄ）と単位時間（Δｔ）とを乗算した値を加算することにより算出する。具体的には、更新部１４は、式（２３）を算出する。
Figure JPOXMLDOC01-appb-M000037
　すなわち、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘｉ（ｔ））と、対象要素の直前時刻（ｔ）における第２変数（ｙｉ（ｔ））とに基づき更新する。
　更新部１４は、Ｓ１０５とＳ１０７との間のループ処理をＮ回実行した場合、処理をＳ１０８に進める。
　続いて、更新部１４は、Ｓ１０８とＳ１１３との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。
　Ｓ１０９において、更新部１４は、Ｎ個の要素のそれぞれの対象時刻（ｔ＋Δｔ）における第１変数（ｘ１（ｔ＋Δｔ）～ｘＮ（ｔ＋Δｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に０－１組合せ最適化問題により予め定められる作用係数と、に基づき更新値（ｚｉ（ｔ＋Δｔ））を算出する。作用係数は、Ｊに含まれる結合係数およびｈに含まれる局所磁場係数である。具体的には、更新部１４は、式（２４）を算出する。
Figure JPOXMLDOC01-appb-M000038
　続いて、Ｓ１１０において、更新部１４は、更新値（ｚｉ（ｔ＋Δｔ））に、係数（ｃ）と－１とを乗算することにより、外力（ｆｉ（ｔ＋Δｔ））を算出する。具体的には、更新部１４は、式（２５）を算出する。
Figure JPOXMLDOC01-appb-M000039
　続いて、Ｓ１１１において、更新部１４は、時間経過に従って増加する関数であるｐ（ｔ＋Δｔ）に基づき定まる値に、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ１（ｔ＋Δｔ））を乗算した時間発展値（ｇｉ（ｔ＋Δｔ））を算出する。具体的には、更新部１４は、式（２６）を算出する。
Figure JPOXMLDOC01-appb-M000040
　続いて、Ｓ１１２において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙｉ（ｔ＋Δ））を、対象要素の直前時刻（ｔ）における第２変数（ｙｉ（ｔ））に、時間発展値（ｇｉ（ｔ＋Δｔ））と外力（ｆｉ（ｔ＋Δｔ））とを加算した値に単位時間（Δｔ）を乗算した値を、加算することにより、算出する。具体的には、更新部１４は、式（２７）を算出する。
Figure JPOXMLDOC01-appb-M000041
　更新部１４は、以上のようなＳ１０８とＳ１１３との間のループ処理をＮ回実行することにより、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第２変数（ｙｉ（ｔ＋Δｔ））を、対象時刻（ｔ＋Δｔ）におけるＮ個の第１変数（ｘ１（ｔ＋Δｔ）～ｘＮ（ｔ＋Δｔ））と、対象要素の直前時刻（ｔ）におけるに第２変数（ｙｉ（ｔ））とに基づき更新する。
　更新部１４は、Ｓ１０８とＳ１１３との間のループ処理をＮ回実行した場合、処理をＳ１１４に進める。
　Ｓ１１４において、更新部１４は、１または複数の制約条件に基づく制約処理を実行する。なお、制約処理については、図５～図７のフローチャートを用いて詳細を後述する。Ｓ１１４の処理を終えると、更新部１４は、処理をＳ１１５に進める。
　Ｓ１１５において、更新部１４は、直前時刻（ｔ）および対象時刻（ｔ＋Δｔ）のそれぞれに単位時間（Δｔ）を加算して、直前時刻（ｔ）および対象時刻（ｔ＋Δｔ）を更新する。Ｓ１１６において、更新部１４は、Ｓ１０５からＳ１１５までの処理を、ｔが終了時刻（Ｔ）を超えるまで繰り返す。そして、更新部１４は、ｔが終了時刻（Ｔ）より大きくなった場合、本フローを終了する。
　そして、出力部１６は、Ｎ個の要素のそれぞれについて、終了時刻（ｔ＝Ｔ）における第１変数（ｘｉ（Ｔ））の符号に応じて、対応するスピンの値を算出する。例えば、出力部１６は、終了時刻（ｔ＝Ｔ）における第１変数（ｘｉ（Ｔ））の符号が負である場合、対応するスピンを－１とし、正である場合、対応するスピンを＋１とする。そして、出力部１６は、算出した複数のスピンの値、または、算出した複数のスピンの値を離散変数に変換した値を組合せ最適化問題の解として出力する。
　以上のＳ１０１～Ｓ１１６の処理を実行することにより、更新部１４は、改良アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））およびＮ個の第２変数（ｙ１（ｔ）～ｙＮ（ｔ））を算出することができる。
　図４は、更新部１４の処理の流れの第２例を示すフローチャートである。更新部１４は、改良アルゴリズムを用いて０－１組合せ最適化問題を解く場合、図３に示す流れに代えて、図４に示す流れで処理を実行してもよい。
　まず、Ｓ２０１、Ｓ２０２およびＳ２０３において、更新部１４は、図３に示す第１例のＳ１０１、Ｓ１０２およびＳ１０３と同一の処理を実行する。
　続いて、更新部１４は、Ｓ２０４とＳ２１６との間のループ処理を、ｔがＴより大きくなるまで繰り返す。１回のループ処理において、更新部１４は、対象時刻（ｔ＋Δｔ）におけるＮ個の第２変数（ｙ１（ｔ＋Δｔ）～ｙＮ（ｔ＋Δｔ））を、直前時刻（ｔ）におけるＮ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））および直前時刻（ｔ）におけるＮ個の第２変数（ｙ１（ｔ）～ｙＮ（ｔ））に基づき算出する。また、１回のループ処理において、更新部１４は、対象時刻（ｔ＋Δｔ）におけるＮ個の第１変数（ｘ１（ｔ＋Δｔ）～ｘＮ（ｔ＋Δｔ））を、直前時刻（ｔ）におけるＮ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））、および、対象時刻（ｔ＋Δｔ）におけるＮ個の第２変数（ｙ１（ｔ＋Δｔ）～ｙＮ（ｔ＋Δｔ））に基づき算出する。
　続いて、更新部１４は、Ｓ２０５とＳ２１０との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ２０５とＳ２１０との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。
　Ｓ２０６において、更新部１４は、Ｎ個の要素のそれぞれの直前時刻（ｔ）における第１変数（ｘ１（ｔ）～ｘＮ（ｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に０－１組合せ最適化問題により予め定められる作用係数と、に基づき更新値（ｚｉ（ｔ））を算出する。具体的には、更新部１４は、式（２８）を算出する。
Figure JPOXMLDOC01-appb-M000042
　続いて、Ｓ２０７において、更新部１４は、更新値（ｚｉ（ｔ））に、係数（ｃ）と－１とを乗算することにより、外力（ｆｉ（ｔ））を算出する。具体的には、更新部１４は、式（２９）を算出する。
Figure JPOXMLDOC01-appb-M000043
　続いて、Ｓ２０８において、更新部１４は、時間経過に従って増加する関数であるｐ（ｔ）に基づき定まる値に、対象要素の直前時刻（ｔ）における第１変数（ｘ１（ｔ））を乗算した時間発展値（ｇｉ（ｔ））を算出する。具体的には、更新部１４は、式（３０）を算出する。
Figure JPOXMLDOC01-appb-M000044
　続いて、Ｓ２０９において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙｉ（ｔ＋Δ））を、対象要素の直前時刻（ｔ）における第２変数（ｙｉ（ｔ））に、時間発展値（ｇｉ（ｔ））と外力（ｆｉ（ｔ））とを加算した値に単位時間（Δｔ）を乗算した値を、加算することにより、算出する。具体的には、更新部１４は、式（３１）を算出する。
Figure JPOXMLDOC01-appb-M000045
　更新部１４は、以上のようなＳ２０５とＳ２１０との間のループ処理をＮ回実行することにより、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第２変数（ｙｉ（ｔ＋Δｔ））を、直前時刻（ｔ）におけるＮ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））と、対象要素の直前時刻（ｔ）におけるに第２変数（ｙｉ（ｔ））とに基づき更新する。
　更新部１４は、Ｓ２０５とＳ２１０との間のループ処理をＮ回実行した場合、処理をＳ２１１に進める。
　Ｓ２１１において、更新部１４は、１または複数の制約条件に基づく制約処理を実行する。なお、制約処理については、図５～図７のフローチャートを用いて詳細を後述する。Ｓ２１１の処理を終えると、更新部１４は、処理をＳ２１２に進める。
　続いて、更新部１４は、Ｓ２１２とＳ２１４との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。
　Ｓ２１３において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘｉ（ｔ））に、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙｉ（ｔ＋Δｔ））と予め定められた定数（Ｄ）と単位時間（Δｔ）とを乗算した値を加算することにより算出する。具体的には、更新部１４は、式（３２）を算出する。
Figure JPOXMLDOC01-appb-M000046
　すなわち、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘｉ（ｔ））と、対象要素の直前時刻（ｔ）における第２変数（ｙｉ（ｔ））とに基づき更新する。
　更新部１４は、Ｓ２１２とＳ２１４との間のループ処理をＮ回実行した場合、処理をＳ２１５に進める。
　Ｓ２１５において、更新部１４は、直前時刻（ｔ）および対象時刻（ｔ＋Δｔ）のそれぞれに単位時間（Δｔ）を加算して、直前時刻（ｔ）および対象時刻（ｔ＋Δｔ）を更新する。Ｓ２１６において、更新部１４は、Ｓ２１２からＳ２１５までの処理を、ｔが終了時刻（Ｔ）を超えるまで繰り返す。そして、更新部１４は、ｔが終了時刻（Ｔ）より大きくなった場合、本フローを終了する。
　そして、出力部１６は、Ｎ個の要素のそれぞれについて、終了時刻（ｔ＝Ｔ）における第１変数（ｘｉ（Ｔ））の符号に応じて、対応するスピンの値を算出する。例えば、出力部１６は、終了時刻（ｔ＝Ｔ）における第１変数（ｘｉ（Ｔ））の符号が負である場合、対応するスピンを－１とし、正である場合、対応するスピンを＋１とする。そして、出力部１６は、算出した複数のスピンの値、または、算出した複数のスピンの値を離散変数に変換した値を組合せ最適化問題の解として出力する。
　以上のＳ２０１～Ｓ２１６の処理を実行することにより、更新部１４は、改良アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ１（ｔ）～ｘＮ（ｔ））およびＮ個の第２変数（ｙ１（ｔ）～ｙＮ（ｔ））を算出することができる。
　（制約処理の第１例）
　図５は、制約処理の第１例の処理の流れを示すフローチャートである。制約条件として、１または複数の制約条件としてＭ個の一次不等式が設定されている場合、更新部１４は、図３のＳ１１４および図４のＳ２１１において、図５に示す流れで処理を実行する。
　まず、更新部１４は、Ｓ３０１とＳ３０８との間のループ処理を、ｍ＝１からｉ＝Ｍまでｍを１ずつインクリメントしながら繰り返す。なお、ｍは、Ｍ個の一次不等式のうちの処理対象となる一次不等式を表すインデックスである。
　Ｓ３０１とＳ３０８とのループ内において、まず、Ｓ３０２において、更新部１４は、制約値（Ｓｍ）を算出する。具体的には、図３のＳ１１４の制約処理において、更新部１４は、式（３３）の演算を実行する。また、図４のＳ２１１の制約処理において、更新部１４は、式（３４）の演算を実行する。
Figure JPOXMLDOC01-appb-M000047
Figure JPOXMLDOC01-appb-M000048
　なお、制御値であるＳｍは、対応する一次不等式における定数項（Ｗｍ）を除く関数に含まれる複数の離散変数のそれぞれに、対応する第１変数（ｘｉ）を代入し得られる値である。
　続いて、Ｓ３０３において、更新部１４は、制約値（Ｓｍ）が対応する一次不等式における定数項（Ｗｍ）以下であるか否かを判断する。すなわち、更新部１４は、複数の第１変数（ｘｉ）が対応する一次不等式を満たすか否かを判断する。
　更新部１４は、制約値（Ｓｍ）が対応する一次不等式における定数項（Ｗｍ）以下である場合、すなわち、複数の第１変数（ｘｉ）が対応する一次不等式を満たす場合（Ｓ３０３のＹｅｓ）、処理をＳ３０８に進める。更新部１４は、制約値（Ｓｍ）が対応する一次不等式における定数項（Ｗｍ）以下ではない場合、すなわち、複数の第１変数（ｘｉ）が対応する一次不等式を満たさない場合（Ｓ３０３のＮｏ）、処理をＳ３０４に進める。
　続いて、更新部１４は、Ｓ３０４とＳ３０７との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ３０４とＳ３０７との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。
　Ｓ３０４とＳ３０７とのループ内において、まず、Ｓ３０５において、更新部１４は、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値（ｒｍ，ｉ）を算出する。具体的には、更新部１４は、式（３５）の演算を実行する。
Figure JPOXMLDOC01-appb-M000049
　続いて、Ｓ３０６において、更新部１４は、対象要素の第２変数（ｙｉ）から、補正値（ｒｍ，ｉ）を減算する。具体的には、更新部１４は、式（３６）の演算を実行する。
Figure JPOXMLDOC01-appb-M000050
　更新部１４は、Ｓ３０４とＳ３０７との間のループ処理をＮ回実行した場合、処理をＳ３０８に進める。
　そして、更新部１４は、Ｓ３０１とＳ３０８との間のループ処理をＭ回実行した場合、すなわち、Ｍ個の一次不等式の全てについて処理を実行した場合、本フローを終了する。
　以上の処理を実行することにより、更新部１４は、Ｍ個の一次不等式のそれぞれについて、複数の第１変数（ｘｉ）が対応する一次不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙｉ）から、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値（ｒｍ，ｉ）を減算することができる。これにより、更新部１４は、複数の第１変数（ｘｉ）が対応する一次不等式を満たしていない場合、複数の第１変数（ｘｉ）が対応する一次不等式を満たす方向に変化するように、複数の第２変数（ｙｉ）を補正することができる。この結果、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。
　（制約処理の第２例）
　図６は、制約処理の第２例の処理の流れを示すフローチャートである。制約条件として、１または複数の制約条件としてＭ個の一次不等式が設定されている場合、更新部１４は、図３のＳ１１４および図４のＳ２１１において、図６に示す流れで処理を実行してもよい。
　まず、更新部１４は、Ｓ４０１とＳ４０９との間のループ処理を、ｍ＝１からｉ＝Ｍまでｍを１ずつインクリメントしながら繰り返す。
　Ｓ４０１とＳ４０９とのループ内において、まず、Ｓ４０２において、更新部１４は、制約値（Ｓｍ）を算出する。Ｓ４０２の処理は、図５のＳ３０２の処理と同一である。
　続いて、Ｓ４０３において、更新部１４は、制約値（Ｓｍ）が対応する一次不等式における定数項（Ｗｍ）以下であるか否かを判断する。Ｓ４０３の処理は、図５のＳ３０３の処理と同一である。更新部１４は、制約値（Ｓｍ）が対応する一次不等式における定数項（Ｗｍ）以下である場合（Ｓ４０３のＹｅｓ）、処理をＳ４０９に進める。更新部１４は、制約値（Ｓｍ）が対応する一次不等式における定数項（Ｗｍ）以下ではない場合（Ｓ４０３のＮｏ）、処理をＳ４０４に進める。
　Ｓ４０４において、更新部１４は、複数の要素により特定される方向における対応する一次不等式の境界に対する傾き（Ｅｍ）を算出する。具体的には、更新部１４は、式（３７）の演算を実行する。
Figure JPOXMLDOC01-appb-M000051
　続いて、更新部１４は、Ｓ４０５とＳ４０８との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ４０５とＳ４０８との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。
　Ｓ４０５とＳ４０８とのループ内において、まず、Ｓ４０６において、更新部１４は、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分、および、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた補正値（ｒｍ，ｉ）を算出する。具体的には、更新部１４は、式（３８）の演算を実行する。
Figure JPOXMLDOC01-appb-M000052
　続いて、Ｓ４０７において、更新部１４は、対象要素の第２変数（ｙｉ）から、補正値（ｒｍ，ｉ）を減算する。Ｓ４０７の処理は、図５のＳ３０６の処理と同一である。
　更新部１４は、Ｓ４０５とＳ４０８との間のループ処理をＮ回実行した場合、処理をＳ４０９に進める。
　そして、更新部１４は、Ｓ４０１とＳ４０９との間のループ処理をＭ回実行した場合、すなわち、Ｍ個の一次不等式の全てについて処理を実行した場合、本フローを終了する。
　以上の処理を実行することにより、更新部１４は、Ｍ個の一次不等式のそれぞれについて、複数の第１変数（ｘｉ）が対応する一次不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙｉ）から、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分、および、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた補正値（ｒｍ，ｉ）を減算することができる。特に、更新部１４は、補正値（ｒｍ，ｉ）に、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた値を加えることにより、単位時間（Δｔ）毎の更新処理において、複数の第１変数（ｘｉ）および複数の第２変数（ｙｉ）が不安定に変動することを抑制することができる。
　これにより、更新部１４は、複数の第１変数（ｘｉ）が対応する一次不等式を満たしていない場合、複数の第１変数（ｘｉ）が対応する一次不等式を満たす方向に安定して変化するように、複数の要素に対応する複数の第２変数（ｙｉ）を補正することができる。この結果、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。
　（制約処理の第３例）
　図７は、制約処理の第３例の処理の流れを示すフローチャートである。１または複数の制約条件のうちの一つに二次不等式が設定されている場合、更新部１４は、二次不等式については、図５のＳ３０２～Ｓ３０７または図６のＳ４０２～Ｓ４０８の処理に代えて、図７に示すＳ５０１～Ｓ５０６の処理を実行する。
　まず、Ｓ５０１において、更新部１４は、制約値（Ｓ）を算出する。具体的には、図３のＳ１１４の制約処理において、更新部１４は、式（３９）の演算を実行する。また、図４のＳ２１１の制約処理において、更新部１４は、式（４０）の演算を実行する。
Figure JPOXMLDOC01-appb-M000053
Figure JPOXMLDOC01-appb-M000054
　なお、制御値であるＳは、二次不等式における定数項（ｑ）を除く関数に含まれる複数の離散変数のそれぞれに、対応する第１変数（ｘｉ）を代入して得られる値である。
　続いて、Ｓ５０２において、更新部１４は、制約値（Ｓ）が二次不等式における定数項（ｑ）以下であるか否かを判断する。すなわち、更新部１４は、複数の第１変数（ｘｉ）が二次不等式を満たすか否かを判断する。
　更新部１４は、制約値（Ｓ）が二次不等式における定数項（ｑ）以下である場合、すなわち、複数の第１変数（ｘｉ）が対応する二次不等式を満たす場合（Ｓ５０２のＹｅｓ）、本フローを終了する。更新部１４は、制約値（Ｓ）が二次不等式における定数項（ｑ）以下ではない場合、すなわち、複数の第１変数（ｘｉ）が対応する二次不等式を満たさない場合（Ｓ５０２のＮｏ）、処理をＳ５０３に進める。
　続いて、更新部１４は、Ｓ５０３とＳ５０６との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ５０３とＳ５０６との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。
　Ｓ５０３とＳ５０６とのループ内において、まず、Ｓ５０４において、更新部１４は、二次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値（ｒｉ）を算出する。具体的には、図３のＳ１１４の制約処理において、更新部１４は、式（４１）の演算を実行する。また、図４のＳ２１１の制約処理において、更新部１４は、式（４２）の演算を実行する。
Figure JPOXMLDOC01-appb-M000055
Figure JPOXMLDOC01-appb-M000056
　続いて、Ｓ５０５において、更新部１４は、対象要素の第２変数（ｙｉ）から、補正値（ｒｉ）を減算する。具体的には、更新部１４は、式（４３）の演算を実行する。
Figure JPOXMLDOC01-appb-M000057
　更新部１４は、Ｓ５０３とＳ５０６との間のループ処理をＮ回実行した場合、本フローを終了する。
　以上の処理を実行することにより、更新部１４は、二次不等式について、複数の第１変数（ｘｉ）が対応する二次不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙｉ）から、二次不等式の境界から複数の要素により特定される位置までの距離の対象要素の成分に応じた補正値（ｒｉ）を減算することができる。これにより、更新部１４は、複数の第１変数（ｘｉ）が対応する二次不等式を満たしていない場合、複数の第１変数（ｘｉ）が対応する二次不等式を満たす方向に変化するように、複数の第２変数（ｙｉ）を補正することができる。この結果、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。
　以上のように、本実施形態に係る求解装置１０は、シミュレーテッド分岐アルゴリズムにより０－１組合せ最適化問題を解く。本実施形態に係る求解装置１０は、単位時間（Δｔ）毎の更新処理において、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数を代入した不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙｉ）から、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値（ｒｉ）を減算する。
　これにより、本実施形態に係る求解装置１０は、不等式により表される１または複数の制約条件の下で、０－１組合せ最適化問題を解く場合であっても、制約条件を設定しない場合と同一の変数の数により、解を得ることができる。従って、本実施形態に係る求解装置１０によれば、スラッグ変数を設けずに、簡易に、不等式により表される１または複数の制約条件の下で０－１組合せ最適化問題を解くことができる。この結果、本実施形態に係る求解装置１０は、変数の途中経過等を記憶するメモリ量を少なくし、計算量および計算時間を少なくすることができる。
　なお、求解装置１０は、一次不等式の制約処理を実行する場合、Ｍ個の一次不等式を生成し、ｉ番目の離散変数に乗算される係数であるωｍ，ｉ、および、一次不等式の定数であるＷｍを呼び出してもよい。
Figure JPOXMLDOC01-appb-M000058
　求解装置１０は、一回の列番号および行番号の指定により、ωｍ，ｉおよびＷｍを同時に読み出すことができ、容易に処理を実行することができる。
　（システム構成）
　図８は、情報処理システム１００の構成を示す図である。改良アルゴリズムは、例えば図８に示すような情報処理システム１００により実行させることが可能である。これにより、情報処理システム１００は、大規模な組合せ最適化問題を、並列処理により高速に解くことができる。
　情報処理システム１００は、管理サーバ１０１と、ネットワーク１０２と、複数の計算サーバ１０３（１０３ａ～１０３ｃ）（情報処理装置）と、複数のケーブル１０４（１０４ａ～１０４ｃ）と、スイッチ１０５を備える。また、端末装置１０６は、情報処理システム１００と通信可能である。管理サーバ１０１、複数の計算サーバ１０３（１０３ａ～１０３ｃ）、端末装置１０６は、ネットワーク１０２を介して互いにデータ通信をする。ネットワーク１０２は、例えば、複数のコンピュータネットワークが相互に接続されたインターネットである。ネットワーク１０２は、通信媒体として有線、無線、または、これらの組合せであってよい。
　また、複数の計算サーバ１０３（１０３ａ～１０３ｃ）は、それぞれケーブル１０４（１０４ａ～１０４ｃ）を介してスイッチ１０５に接続される。複数のケーブル１０４（１０４ａ～１０４ｃ）およびスイッチ１０５は、計算サーバ１０３と計算サーバ１０３と間のインターコネクトを形成する。複数の計算サーバ１０３（１０３ａ～１０３ｃ）のそれぞれは、インターコネクトを介して互いにデータ通信をすることも可能である。スイッチ１０５は、例えば、Ｉｎｆｉｎｉｂａｎｄのスイッチでる。ケーブル１０４ａ～１０４ｃは、例えば、Ｉｎｆｉｎｉｂａｎｄのケーブルである。ただし、スイッチ１０５およびケーブル１０４は、有線ＬＡＮのスイッチ／ケーブルであってもよい。ケーブル１０４およびスイッチ１０５で使われる通信規格および通信プロトコルについては、特に問わない。端末装置１０６は、例えば、ノートＰＣ、デスクトップＰＣ、スマートフォン、タブレットまたは車載端末装置である。
　改良アルゴリズムを用いた組合せ最適化問題の求解処理は、並列的および分散的に実行可能である。従って、複数の計算サーバ１０３（１０３ａ～１０３ｃ）のそれぞれおよび／または計算サーバ１０３（１０３ａ～１０３ｃ）のプロセッサは、改良アルゴリズムを用いた組合せ最適化問題の求解処理における一部の計算処理の一部のステップを分担して実行してもよいし、異なる変数に対する同一の計算処理を並列的に実行してもよい。この場合、管理サーバ１０１は、例えば、ユーザによって入力された組合せ最適化問題を各計算サーバ１０３に処理可能な形式に変換し、各計算サーバ１０３に実行させる。そして、管理サーバ１０１は、各計算サーバ１０３から計算結果を取得し、集約した計算結果を組合せ最適化問題の解に変換する。
　図８には、３台の計算サーバ１０３（１０３ａ～１０３ｃ）が示されている。しかし、情報処理システム１００に含まれる計算サーバ１０３の台数は、３台に限定されない。例えば、情報処理システム１００は、１台であっても、２台であっても、４台以上であってもよい。また、情報処理システム１００は、含んでいる複数の計算サーバ１０３のうちの一部を用いて、組合せ最適化問題の求解の処理を実行してもよい。計算サーバ１０３は、どのような種類の情報処理装置であってもよい。例えば、計算サーバ１０３は、データセンターに設置されたサーバであってもよいし、オフィスに設置されたデスクトップＰＣであってもよい。また、計算サーバ１０３は、異なるローケーションに設置された複数の種類のコンピュータであってもよい。例えば、計算サーバ１０３は、汎用的なコンピュータであってもよいし、専用の電子回路または、これらの組合せであってもよい。
　図９は、管理サーバ１０１の構成を示す図である。管理サーバ１０１は、例えば、中央演算処理装置（ＣＰＵ）とメモリとを含むコンピュータである。管理サーバ１０１は、プロセッサ１１０と、記憶部１１４と、通信回路１１５と、入力回路１１６と、出力回路１１７とを備える。プロセッサ１１０、記憶部１１４、通信回路１１５、入力回路１１６、出力回路１１７は、互いにバス１２０を介して接続される。プロセッサ１１０は、内部の機能構成として、管理部１１１と、変換部１１２と、制御部１１３を含む。
　プロセッサ１１０は、演算を実行し、管理サーバ１０１の制御を行う電子回路である。プロセッサ１１０は、例えば、ＣＰＵ、マイクロプロセッサ、ＡＳＩＣ、ＦＰＧＡ、ＰＬＤまたはこれらの組合せであってよい。管理部１１１は、ユーザの端末装置１０６を介して管理サーバ１０１の操作を行うためのインタフェースを提供する。管理部１１１が提供するインタフェースは、例えば、ＡＰＩ、ＣＬＩまたはウェブページ等である。例えば、ユーザは、管理部１１１を介して組合せ最適化問題の情報の入力を行ったり、計算された組合せ最適化問題の解の閲覧および／またはダウンロードを行ったりする。変換部１１２は、組合せ最適化問題に関するパラメータを入力して、入力したパラメータを各計算サーバ１０３が処理可能な形式に変換する。制御部１１３は、各計算サーバ１０３に制御指令を送信する。制御部１１３が各計算サーバ１０３から計算結果を取得した後、変換部１１２は、複数の計算結果を集約し、組合せ最適化問題の解に変換し、組合せ最適化問題の解を出力する。
　記憶部１１４は、管理サーバ１０１のプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。プログラムは、ＯＳとアプリケーションの両方を含む。記憶部１１４は、揮発性メモリ、不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリは、例えば、ＤＲＡＭまたはＳＲＡＭ等である。不揮発性メモリは、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭまたはＭＲＡＭ等である。また、記憶部１１４は、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置等であってもよい。
　通信回路１１５は、ネットワーク１０２に接続された各装置との間でデータの送受信を行う。通信回路１１５は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ　Ｉｎｔｅｒｆａｃｅ　Ｃａｒｄ）である。ただし、通信回路１１５は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。入力回路１１６は、管理サーバ１０１へのデータ入力を実現する。入力回路１１６は、外部ポートとして、例えば、ＵＳＢ、ＰＣＩ－Ｅｘｐｒｅｓｓなどを含む。操作装置１１８は、入力回路１１６に接続される。操作装置１１８は、ユーザが管理サーバ１０１に情報を入力するための装置である。操作装置１１８は、例えば、キーボード、マウス、タッチパネル、音声認識装置などであるが、これらに限られない。出力回路１１７は、管理サーバ１０１からのデータ出力を実現する。出力回路１１７は、外部ポートとしてＨＤＭＩ（登録商標）またはＤｉｓｐｌａｙＰｏｒｔ等を含む。例えば、表示装置１１９は、出力回路１１７に接続される。表示装置１１９は、例えば、ＬＣＤ（液晶ディスプレイ）、有機ＥＬ（有機エレクトロルミネッセンス）ディスプレイまたはプロジェクタ等であるが、これらに限られない。
　管理者は、操作装置１１８および表示装置１１９を使って、管理サーバ１０１のメンテナンスを行うことができる。なお、操作装置１１８および表示装置１１９は、管理サーバ１０１に組み込まれたものであってもよい。また、操作装置１１８および表示装置１１９は、管理サーバ１０１に接続されていなくてもよい。例えば、管理者は、ネットワーク１０２と通信可能な端末装置を用いて管理サーバ１０１のメンテナンスを行ってもよい。
　図１０は、管理サーバ１０１の記憶部１１４に記憶されるデータを示す図である。記憶部１１４は、例えば、問題データ１１４Ａと、計算データ１１４Ｂと、管理プログラム１１４Ｃと、変換プログラム１１４Ｄと、制御プログラム１１４Ｅとを記憶する。例えば、問題データ１１４Ａは、組合せ最適化問題のデータを含む。例えば、計算データ１１４Ｂは、各計算サーバ１０３から収集された計算結果を含む。例えば、管理プログラム１１４Ｃは、管理部１１１の機能を実現するプログラムである。例えば、変換プログラム１１４Ｄは、変換部１１２の機能を実現するプログラムである。例えば、制御プログラム１１４Ｅは、制御部１１３の機能を実現するプログラムである。
　図１１は、計算サーバ１０３ａの構成を示す図である。他の計算サーバ１０３は、計算サーバ１０３ａと同様の構成であってもよいし、計算サーバ１０３ａと異なる構成であってもよい。
　計算サーバ１０３ａは、例えば、通信回路１３１と、共有メモリ１３２と、プロセッサ１３３ａ～１３３ｄと、ストレージ１３４と、ホストバスアダプタ１３５とを備える。通信回路１３１、共有メモリ１３２、プロセッサ１３３ａ～１３３ｄ、ストレージ１３４、ホストバスアダプタ１３５は、バス１３６を介して互いに接続される。
　通信回路１３１は、ネットワーク１０２に接続された各装置との間でデータの送受信を行う。通信回路１３１は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ　Ｉｎｔｅｒｆａｃｅ　Ｃａｒｄ）である。ただし、通信回路１３１は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。共有メモリ１３２は、プロセッサ１３３ａ～１３３ｄからアクセス可能なメモリである。共有メモリ１３２は、例えば、ＤＲＡＭおよびＳＲＡＭ等の揮発性メモリである。共有メモリ１３２は、不揮発性メモリ等の他の種類のメモリを含んでもよい。プロセッサ１３３ａ～１３３ｄは、共有メモリ１３２を介してデータを共有する。なお、共有メモリ１３２は、計算サーバ１０３ａの全てメモリにより構成されていなくてもよい。例えば、計算サーバ１０３ａの一部のメモリは、プロセッサ１３３ａ～１３３ｄのうちのいずれかからのみからアクセスできるローカルメモリであってもよい。
　プロセッサ１３３ａ～１３３ｄは、計算処理を実行する電子回路である。プロセッサ１３３ａ～１３３ｄは、例えば、ＣＰＵ、ＧＰＵ、ＦＰＧＡ、ＡＳＩＣのいずれであってもよいし、これらの組合せであってもよい。また、プロセッサ１３３ａ～１３３ｄは、ＣＰＵコアまたはＣＰＵスレッドであってもよい。プロセッサ１３３ａ～１３３ｄがＣＰＵである場合、計算サーバ１０３ａが備えるソケット数については、特に問わない。また、プロセッサ１３３ａ～１３３ｄは、ＰＣＩ　ｅｘｐｒｅｓｓなどのバスを介して計算サーバ１０３ａのその他の構成要素に接続されていてもよい。
　図１１の例では、計算サーバ１０３ａは、４つのプロセッサ１３３ａ～１３３ｄを備える。しかし、１台の計算サーバ１０３ａに含まれるプロセッサ数は、４個に限られない。
　ストレージ１３４は、計算サーバ１０３ａのプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、ＯＳとアプリケーションの両方を含むものとする。ストレージ１３４は、揮発性メモリおよび不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリは、たとえば、ＤＲＡＭまたはＳＲＡＭである。不揮発性メモリは、例えば、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭまたはＭＲＡＭ等である。また、ストレージ１３４は、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置を含んでもよい。
　ホストバスアダプタ１３５は、他の計算サーバ１０３との間のデータ通信を実現する。ホストバスアダプタ１３５は、ケーブル１０４ａを介してスイッチ１０５に接続されている。ホストバスアダプタ１３５は、例えば、ＨＣＡ（Ｈｏｓｔ　Ｃｈａｎｎｅｌ　Ａｄａｐｔｏｒ）である。ホストバスアダプタ１３５、ケーブル１０４ａ、スイッチ１０５で高スループットを実現可能なインターコネクトを形成することにより、並列的な計算処理の速度を向上させることができる。
　図１２は、ストレージ１３４に記憶されるデータを示す図である。ストレージ１３４は、例えば、計算データ１３４Ａと、計算プログラム１３４Ｂと、制御プログラム１３４Ｃとを記憶する。計算データ１３４Ａは、計算サーバ１０３ａの計算途中のデータまたは計算結果を含む。計算データ１３４Ａの少なくとも一部は、共有メモリ１３２、プロセッサのキャッシュまたはプロセッサのレジスタ等の、異なる記憶階層に記憶されてもよい。計算プログラム１３４Ｂは、各プロセッサ１３３における計算処理および、共有メモリ１３２およびストレージ１３４へのデータの保存処理を実現するプログラムである。制御プログラム１３４Ｃは、管理サーバ１０１の制御部１１３から送信された指令に基づき、計算サーバ１０３ａを制御し、計算サーバ１０３ａの計算結果を管理サーバ１０１に送信するプログラムである。
　このような計算サーバ１０３ａは、プロセッサ１３３ａ～１３３ｄが組合せ最適化問題を解くためのプログラムを実行する。このプログラムは、計算サーバ１０３ａを、入力部１２、更新部１４および出力部１６として機能させる。
　このような構成の情報処理システム１００は、ＰＣクラスタとして利用することが可能である。ＰＣクラスタは、複数台のコンピュータを接続し、１台のコンピュータでは得られない計算性能を実現するシステムである。情報処理システム１００は、例えばＭＰＩ（Ｍｅｓｓａｇｅ　Ｐａｓｓｉｎｇ　Ｉｎｔｅｒｆａｃｅ）を使うことにより、複数の計算サーバ１０３にメモリが分散して配置されている構成であっても、並列的な計算を実行することが可能である。
　なお、情報処理システム１００は、高速リンクで接続された複数のＧＰＵであってもよい。この場合、複数のＧＰＵのそれぞれは、計算サーバ１０３と同様の処理を実行する。
　本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組合せにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組合せてもよい。
１０　求解装置
１２　入力部
１４　更新部
１６　出力部

Claims (16)

  1. 　０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く求解装置であって、
    　第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
    　前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する出力部と、
    　を備え、
    　前記複数の要素は、前記複数の離散変数に対応し、
    　前記第１変数および前記第２変数のそれぞれは、実数により表され、
    　前記単位時間毎の更新処理において、前記更新部は、
    　前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、
    　前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新し、
    　前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
    　求解装置。
  2. 　前記出力部は、
    　前記終了時刻における前記複数の要素のそれぞれについて、前記第１変数を予め設定されたしきい値により２値化した離散変数の値を算出し、
    　算出した前記複数の離散変数の値を前記０－１組合せ最適化問題の解として出力する
    　請求項１に記載の求解装置。
  3. 　前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
    　前記第１変数を更新した後に前記第２変数を更新し、
    　前記第２変数を更新した後に前記第２変数から前記補正値を減算する
    　請求項２に記載の求解装置。
  4. 　前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
    　前記第２変数を更新した後に前記第１変数を更新し、
    　前記第２変数を更新した後、且つ、前記第１変数を更新する前に、前記第２変数から前記補正値を減算する
    　請求項２に記載の求解装置。
  5. 　前記単位時間毎の更新処理において、前記更新部は、
    　前記複数の要素のそれぞれについて、対象時刻における前記第１変数を、前記対象時刻より単位時間前の直前時刻における前記第１変数に、前記第２変数と予め定められた定数と前記単位時間とを乗算した値を加算することにより、算出する
    　請求項２に記載の求解装置。
  6. 　前記単位時間毎の更新処理において、前記更新部は、
    　前記複数の要素のそれぞれについて、
    　前記複数の要素のそれぞれの前記第１変数と、対象要素と前記複数の要素のそれぞれとの組毎に前記０－１組合せ最適化問題により定められる作用係数とに基づき、外力を算出し、
    　時間経過に従って増加する関数に基づき定まる値と、前記対象要素の前記第１変数とを乗算した時間発展値を算出し、
    　前記直前時刻における前記第２変数に、前記時間発展値と前記外力とを加算した値に前記単位時間を乗算した値を、加算することにより、前記対象時刻における前記第２変数を算出する
    　請求項５に記載の求解装置。
  7. 　前記０－１組合せ最適化問題は、Ｎ個の離散変数を含み、
    　前記更新部は、前記Ｎ個の離散変数のうちのｉ番目の離散変数に対応するｉ番目の要素の、前記対象時刻における前記第１変数を式（１０１）または式（１０２）により算出し、
    Figure JPOXMLDOC01-appb-M000001
     　Ｎは、２以上の整数であり、
    　ｉは、１からＮまでの整数であり、
    　Ｄは、前記定数であり、
    　Δｔは、前記単位時間であり、
    　ｔは、前記直前時刻であり、
    　ｔ＋Δｔは、前記対象時刻であり、
    　ｘｉ（ｔ）は、前記ｉ番目の要素の前記直前時刻における前記第１変数であり、
    　ｙｉ（ｔ）は、前記ｉ番目の要素の前記直前時刻における前記第２変数であり、
    　ｘｉ（ｔ＋Δｔ）は、前記ｉ番目の要素の前記対象時刻における前記第１変数であり、
    　ｙｉ（ｔ＋Δｔ）は、前記ｉ番目の要素の前記対象時刻における前記第２変数である
    　請求項６に記載の求解装置。
  8. 　前記０－１組合せ最適化問題は、０－１二次計画問題であり、
    　前記更新部は、前記ｉ番目の要素の前記対象時刻における前記第２変数を式（１０３）または式（１０４）により算出し、
    Figure JPOXMLDOC01-appb-M000002
     　ｇｉ（ｔ）およびｇｉ（ｔ＋Δｔ）は、前記時間発展値であり、
    　ｆｉ（ｔ＋Δｔ）は、式（１０５）により表され、ｆｉ（ｔ）は、式（１０６）により表され、
    Figure JPOXMLDOC01-appb-M000003
     　ｚｉ（ｔ＋Δｔ）は、式（１０７）により表され、
    　ｚｉ（ｔ）は、式（１０８）により表され、
    Figure JPOXMLDOC01-appb-M000004
     　ｊは、１からＮまでの整数であり、
    　ｈｉは、Ｎ個の局所磁場係数を含む予め定められた配列に含まれる、ｉ番目の局所磁場係数であり、
    　Ｊｉ，ｊは、Ｎ×Ｎ個の結合係数を含む予め定められた行列に含まれる、ｉ行、ｊ列の結合係数であり、
    　ｃは、予め定められた実数であり、
    　ｘｊ（ｔ）は、前記Ｎ個の離散変数のうちのｊ番目の離散変数に対応するｊ番目の要素の前記直前時刻における前記第１変数であり、
    　ｘｊ（ｔ＋Δｔ）は、前記ｊ番目の要素の前記対象時刻における前記第１変数であり、
    　α（ｔ）は、ｔを変数とする予め定められた関数である
    　請求項７に記載の求解装置。
  9. 　ｇｉ（ｔ）は、式（１０９）により表され、
    　ｇｉ（ｔ＋Δｔ）は、式（１１０）により表され、
    Figure JPOXMLDOC01-appb-M000005
     　ｐ（ｔ）は、ｔを変数とする予め定められた関数であり、ｔに従って増加し、ｔが前記初期時刻において０となり、ｔが前記終了時刻において１となり、
    　Ｋは、予め定められた定数である
    　請求項８に記載の求解装置。
  10. 　前記１または複数の制約条件のうちのｍ番目の制約条件は、式（１１１）により表され、
    Figure JPOXMLDOC01-appb-M000006
     　ｍは、１以上の整数であり、
    　ｘｉは、ｉ番目の離散変数であり、
    　ωｍ，ｉは、ｉ番目の離散変数に乗算される係数であり、
    　Ｗｍは、予め定められた定数である
    　請求項８または９に記載の求解装置。
  11. 　前記更新部は、前記ｉ番目の要素の前記第２変数に対する前記補正値として、式（１１２）により表されるｒｍ，ｉを算出し、
    Figure JPOXMLDOC01-appb-M000007
     　Ａｍは、前記ｍ番目の制約条件に対して予め定められた係数であり、
    　ｋは、予め定められた１以上の整数であり、
    　Ｓｍは、式（１１３）または式（１１４）により表される
    Figure JPOXMLDOC01-appb-M000008
     　請求項１０に記載の求解装置。
  12. 　前記更新部は、前記ｉ番目の要素の前記第２変数に対する前記補正値として、式（１１５）により表されるｒｍ，ｉを算出し、
    Figure JPOXMLDOC01-appb-M000009
     　Ａｍは、前記ｍ番目の制約条件に対して予め定められた係数であり、
    　ｋは、予め定められた１以上の整数であり、
    　Ｓｍは、式（１１６）または式（１１７）により表され、
    Figure JPOXMLDOC01-appb-M000010
     　Ｂｍは、０以上１以下の値であり、前記ｍ番目の制約条件に対して予め定められた係数であり、
    　Ｅｍは、式（１１８）により表される
    Figure JPOXMLDOC01-appb-M000011
     　請求項１０に記載の求解装置。
  13. 　前記１または複数の制約条件のうちの一つの制約条件は、式（１１９）により表され、
    Figure JPOXMLDOC01-appb-M000012
     　ｘｉは、ｉ番目の離散変数であり、
    　ｘｊは、ｊ番目の離散変数であり、
    　Ｑｉ，ｊは、半正定値行列に含まれるｉ行、ｊ列の値であり、
    　ｑは、予め定められた定数である
    　請求項８または９に記載の求解装置。
  14. 　前記更新部は、前記ｉ番目の要素の前記第２変数に対する前記補正値として、式（１２０）または（１２１）により表されるｒｉを算出し、
    Figure JPOXMLDOC01-appb-M000013
     　Ａは、対応する制約条件に対して予め定められた係数であり、
    　ｋは、予め定められた１以上の整数であり、
    　Ｓは、式（１２０）の場合は式（１２２）により表され、式（１２１）の場合には式（１２３）により表される
    Figure JPOXMLDOC01-appb-M000014
     　請求項１３に記載の求解装置。
  15. 　情報処理装置により、０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く求解方法であって、
    　前記情報処理装置により、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新ステップと、
    　前記情報処理装置により、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する出力ステップと、
    　を含み、
    　前記複数の要素は、前記複数の離散変数に対応し、
    　前記第１変数および前記第２変数のそれぞれは、実数により表され、
    　前記単位時間毎の更新ステップにおいて、前記情報処理装置は、
    　前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、
    　前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新し、
    　前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
    　求解方法。
  16. 　情報処理装置を、０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く求解装置として機能させるためのプログラムであって、
    　前記情報処理装置を、
    　第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
    　前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する出力部と、
    　して機能させ、
    　前記複数の要素は、前記複数の離散変数に対応し、
    　前記第１変数および前記第２変数のそれぞれは、実数により表され、
    　前記単位時間毎の更新処理において、前記更新部は、
    　前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、
    　前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新し、
    　前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
    　プログラム。
PCT/JP2022/017778 2021-05-24 2022-04-14 求解装置、求解方法およびプログラム WO2022249785A1 (ja)

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
JP2021-087125 2021-05-24
JP2021087125A JP2022180174A (ja) 2021-05-24 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)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020196862A1 (ja) * 2019-03-28 2020-10-01 株式会社 東芝 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020196862A1 (ja) * 2019-03-28 2020-10-01 株式会社 東芝 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
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
CA3219868A1 (en) 2022-12-01
US20240095300A1 (en) 2024-03-21
JP2022180174A (ja) 2022-12-06

Similar Documents

Publication Publication Date Title
US20210224447A1 (en) Grouping of pauli strings using entangled measurements
US20220012387A1 (en) Information processing device, information processing system, information processing method, and storage medium
JP7297540B2 (ja) 情報処理装置、ｐｕｂｏソルバ、情報処理方法およびプログラム
CN113379058B (zh) 量子模拟方法及装置、电子设备及存储介质
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
Leleu et al. Scaling advantage of nonrelaxational dynamics for high-performance combinatorial optimization
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