Calculation device, calculation method, and computer program product

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is based upon and claims the benefit of priority from Japanese Patent Application No. 2021-036688, filed on Mar. 8, 2021; the entire contents of which are incorporated herein by reference.

FIELD

Embodiments described herein relate generally to a calculation device, a calculation method, and a computer program product.

BACKGROUND

Combinatorial optimization problems are problems of selecting a combination most suitable for a purpose from among a plurality of combinations. Combinatorial optimization problems mathematically result in problems of maximizing a function, called an “objective function”, having a plurality of discrete variables, or problems of minimizing the function. While combinatorial optimization problems are universal problems in various fields such as finance, logistics, transportation, designing, manufacturing, and life science, optimal solutions are not always successfully calculated because of “combinatorial explosion” in which the number of combinations increases exponentially with the problem size. Moreover, it is often difficult to obtain even an approximate solution close to the optimal solution.

There has been a demand for developing technologies for accurately calculating a solution to a combinatorial optimization problem in order to solve problems in various fields and promote social innovation and the progress in science technologies.

An object of an embodiment is to provide a calculation device, a calculation method, and a computer program product to accurately calculate a solution to a combinatorial optimization problem.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram illustrating a bifurcation phenomenon of a simulated bifurcation algorithm;

FIG. 2 is a diagram illustrating a bifurcation phenomenon of an improved algorithm;

FIG. 3 is a diagram illustrating a mean value of the numbers of cuts obtained when a maximum cut problem was solved;

FIG. 4 is a diagram illustrating a maximum value of the numbers of cuts obtained when the maximum cut problem was solved;

FIG. 5 is a diagram illustrating a success probability obtained when the maximum cut problem was solved;

FIG. 6 is a functional configuration diagram of a calculation device according to the present embodiment;

FIG. 7 is a flowchart illustrating a first example of a process of an updating module;

FIG. 8 is a flowchart illustrating a second example of a process of the updating module;

FIG. 9 is a block configuration diagram of the updating module;

FIG. 10 is a configuration diagram of an information processing system;

FIG. 11 is a configuration diagram of a management server;

FIG. 12 is a diagram illustrating data stored in a storage unit;

FIG. 13 is a configuration diagram illustrating a calculation server; and

FIG. 14 is a diagram illustrating data stored in a storage.

DETAILED DESCRIPTION

According to an embodiment, a calculation device configured to solve a combinatorial optimization problem includes a memory and one or more processors coupled to the memory. The one or more processors are configured to: alternately update, for a plurality of elements each associated with a first variable and a second variable, the first variable and the second variable, sequentially for unit times from an initial time to an end time; and output a solution to the combinatorial optimization problem based on the first variables of the plurality of elements at the end time. The plurality of elements correspond to a plurality of discrete variables of the combinatorial optimization problem. The first variables and the second variables are each represented by a real number. The one or more processors are configured to: in an updating process for each of the unit times, for each of the plurality of elements, update the first variable based on the second variable; when the first variable is smaller than a predetermined first value, change the first variable to the first value and change the second variable to a predetermined third value; when the first variable is greater than a predetermined second value greater than the first value, change the first variable to the second value and change the second variable to the third value; and add an acceleration value calculated by a predetermined computation to the second variable.

Embodiments will be described below.

Combinatorial Optimization Problem

An Ising machine is an example of devices used for solving the Ising problem. The Ising machine calculates energy of the ground state of the Ising model. So far, the Ising model has been often used mainly as a model of ferromagnetic or phase transition phenomena. However, in recent years, the Ising model is increasingly used as a model for solving a combinatorial optimization problem. Equation (1) below shows energy of the Ising model.

$\begin{matrix} E_{Ising} = - \sum_{i = 1}^{N} h_{i} s_{i} - \frac{1}{2} \sum_{i = 1}^{N} \sum_{j = 1}^{N} J_{i, j} s_{i} s_{j} & (1) \end{matrix}$

Here, s_iand s_jare spins. The spin is a binary variable that takes a value of either +1 or −1. s_iis the i^thspin. s_jis the j^thspin. i and j are any integer of 1 through N, both inclusive. N denotes the number of spins and is an integer equal to or greater than 2. h_idenotes a local magnetic field acting on the i^thspin. J is a matrix of a coupling coefficient representing a force acting between two spins. J is a real symmetric matrix in which diagonal components are zero. J_ijdenotes an element at the i^throw and j^thcolumn of J. In other words, J_ijis the coupling coefficient representing a force acting between the i^thspin and the j^thspin. The Ising model in Equation (1) is a quadratic equation of spins. The Ising model may be an extended model including a term of 3 or more orders for spins (the Ising model having multibody interaction).

The Ising machine sets the energy E_Isingwritten by Equation (1) as an objective function and calculates a solution that can minimize the energy E_Ising. The solution of the Ising model (s₁, s₂, . . . , s_N) in which the energy E_Isinghas the minimum value is called optimal solution. However, the solution of the Ising model is not necessarily an optimal solution but may be an approximate solution in which the energy E_Isingis close to the minimum value. In other words, the Ising problem may be a problem to calculate not only an optimal solution but also an approximate solution.

A combinatorial optimization problem with an objective function which is a quadratic function with a discrete variable (bit) taking either 0 or 1 is called a quadratic unconstrained binary optimization (QUBO) problem. The discrete variable (bit) is converted into s_iby using a computation of (1+s_i)/2. That is, it can be said that the QUBO problem is equivalent to the Ising problem given by Equation (1). Therefore, the 0-1 combinatorial optimization problem can be converted into the Ising problem to be solved by the Ising machine.

The Ising machine is implemented by hardware, for example, by a quantum annealer, a coherent Ising machine, and a quantum bifurcation machine. The quantum annealer implements quantum annealing using superconducting circuitry. The coherent Ising machine uses the oscillation phenomenon of a network formed with an optical parametric oscillator. The quantum bifurcation machine uses a quantum-mechanical bifurcation phenomenon in a network of a parametric oscillator having Kerr effect. While the Ising machine implemented by such hardware may achieve significant reduction of computation time, scale increase and stable operation are difficult.

A solution to the Ising problem may be calculated using widespread digital computers. Digital computers can be increased in scale and stably operated, compared with quantum annealers, coherent Ising machines, and quantum bifurcation machines. Simulated annealing (SA) is an example of algorithms for solving the Ising problem using a digital computer. However, since simulated annealing is a sequentially updating algorithm in which individual variables are sequentially updated, it is difficult to accelerate a calculation process through parallelization.

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 develops a simulated bifurcation algorithm as an algorithm for solving a combinatorial optimization problem. The simulated bifurcation algorithm can solve a large-scale combinatorial optimization problem with a digital computer using the Ising model. The simulated bifurcation algorithm also can solve a large-scale combinatorial optimization problem with electronic circuitry such as a central processing unit (CPU), a microprocessor, a graphics processing unit (GPU), a field-programmable gate array (FPGA), an application specific integrated circuit (ASIC), or a combination circuit thereof.

Simulated Bifurcation Algorithm The simulated bifurcation algorithm uses two variables x_iand y_ieach corresponding to N elements. The variable x_imay be called first variable, and the variable y_imay be called second variable. In the simulated bifurcation algorithm, each of N elements represents a virtual particle. N elements correspond to N spins in the Ising problem. Therefore, N elements correspond to N discrete variables (bits) of a combinatorial optimization problem. Both the variable x_iand the variable y_iare continuous variables represented by real numbers. The variable x_idenotes the position of the i^thparticle of N particles. The variable y_idenotes the momentum of the i^thparticle. N denotes the number of spins included in the Ising model and is an integer equal to or greater than 2, and i denotes any integer of 1 through N, both inclusive, and denotes an index that identifies each of N elements.

The simulated bifurcation algorithm numerically solves simultaneous ordinary differential equations in Equation (2) below for N variables x_iand N variables y_i.

$\begin{matrix} \frac{{dx}_{i}}{dt} = \frac{\partial H}{\partial y_{i}} = {Dy}_{i} \frac{{dy}_{i}}{dt} = - \frac{\partial H}{\partial x_{i}} = {- D + p (t) - {Kx}_{i}^{2}} x_{i} + f_{i} & (2) \end{matrix}$

H is the Hamiltonian in Equation (3) below.

$\begin{matrix} H = \sum_{i = 1}^{N} [\frac{D}{2} (x_{i}^{2} + y_{i}^{2}) - \frac{p (t)}{2} x_{i}^{2} + \frac{K}{4} x_{i}^{4} - c (h_{i} x_{i} α (t) + \frac{1}{2} \sum_{j = 1}^{N} J_{i, j} x_{i} x_{j})] & (3) \end{matrix}$

The coefficient D is a predetermined constant and corresponds to detuning. The coefficient p(t) corresponds to pumping amplitude and its value monotonously increases according to the number of times of updating in calculation of the simulated bifurcation algorithm. The variable t represents time. The initial value of the coefficient p(t) may be set to 0. The coefficient K corresponds to positive Kerr coefficient.

f_iis external force and given by Equation (4) below.

$\begin{matrix} f_{i} = - {cz}_{i} z_{i} = - h_{i} α (t) - \sum_{j = 1}^{N} J_{i, j} x_{j} & (4) \end{matrix}$

In Equation (4), z_iis given by a partial derivative of the expression inside the parentheses in Equation (3) with respect to the variable x_i. The expression inside the parentheses in Equation (3) corresponds to the energy E_Isingof the Ising model.

The coefficient c may be, for example, a constant determined in advance before calculation is performed. For example, c may be set to a value close to the reciprocal of the maximum eigenvalue of the matrix J. For example, c may be 0.5D√(N/2n). Here, n is the number of edges of a graph for a combinatorial optimization problem. Furthermore, α(t) is a coefficient increasing with p(t). For example, α(t) may be √(p(t)).

The simulated bifurcation algorithm can solve a combinatorial optimization problem to minimize an objective function including a discrete variable of 3 or more orders. The combinatorial optimization problem to minimize an objective function including a discrete variable of 3 or more orders is called a higher order binary optimization (HOBO) problem. To solve a HOBO problem, the simulated bifurcation algorithm uses the Ising model extended to higher order.

The energy in the Ising model extended to higher order is given by Equation (5) below.

$\begin{matrix} E_{HOBO} = \sum_{i = 1}^{N} J_{i}^{(1)} s_{i} + \frac{1}{2} \sum_{i = 1}^{N} \sum_{j = 1}^{N} J_{i, j}^{(2)} s_{i} s_{j} + \frac{1}{3!} \sum_{i = 1}^{N} \sum_{j = 1}^{N} \sum_{k = 1}^{N} J_{i, j, k}^{(3)} s_{i} s_{j} s_{k} + \dots & (5) \end{matrix}$

Here, J⁽ⁿ⁾is an n^th-rank tensor. J⁽ⁿ⁾is generalization of the local magnetic field h_iand the matrix J of the coupling coefficient in Equation (1). For example, J^(l)corresponds to h_i. In J⁽ⁿ⁾, when a plurality of subscripts have the same value, the value of an element is 0. Although Equation (5) denotes the terms up to order 3, order 4 and higher can be defined similarly. Equation (5) corresponds to the energy of the Ising model including a multibody interaction.

When the simulated bifurcation algorithm solves a HOBO problem, the Hamiltonian H in Equation (3) can be replaced by Equation (6).

$\begin{matrix} H = \sum_{i = J}^{N} [\frac{D}{2} (x_{i}^{2} + y_{i}^{2}) - \frac{p (t)}{2} x_{i}^{2} + \frac{K}{4} x_{i}^{4} + c (J_{i}^{(1)} x_{i} α (t) + \frac{1}{2} \sum_{j = 1}^{N} J_{i, j}^{(2)} x_{i} x_{j} + \frac{1}{3!} \sum_{j = 1}^{N} \sum_{k = 1}^{N} J_{i, j, k}^{(3)} x_{i} x_{j} x_{k} + \dots)] H_{Ising} = \sum_{i = 1}^{N} [J_{i}^{(1)} x_{i} α (t) + \frac{1}{2} \sum_{j = 1}^{N} J_{i, j}^{(2)} x_{i} x_{j} + \frac{1}{3!} \sum_{j = 1}^{N} \sum_{k = 1}^{N} J_{i, j, k}^{(3)} x_{i} x_{j} x_{k} + \dots] & (6) \end{matrix}$

When the simulated bifurcation algorithm solves a HOBO problem, the external force f_iin Equation (4) is replaced by Equation (7).

$\begin{matrix} f_{i} = - {cz}_{i} z_{i} = \frac{\partial H_{Ising}}{\partial x_{i}} = J_{i}^{(1)} α (t) + \sum_{j = 1}^{N} J_{i, j}^{(2)} x_{j} + \sum_{j = 1}^{N} \sum_{k = 1}^{N} J_{i, j, k}^{(3)} x_{j} x_{k} + \dots & (7) \end{matrix}$

In Equation (7), z_iis given by a partial derivative of the expression inside the parentheses in Equation (6) with respect to the variable x_i. The expression inside the parentheses in Equation (6) corresponds to the energy E_Isingof the Ising model extended to higher order.

The simulated bifurcation algorithm then calculates the value of the spin s_ibased on the sign of the variable x_iafter the value of p(t) is increased from an initial value (for example, 0) to a prescribed value. The simulated bifurcation algorithm calculates the value of the spin s_iusing the signum function, for example, in which when x_i>0, sgn(x_i)=1, and when x_i<0, sgn(x_i)=−1.

Computation of Simulated Bifurcation Algorithm

When a QUBO problem is solved, the simulated bifurcation algorithm solves differential equations given by Equation (2), Equation (3), and Equation (4) using the symplectic Euler method. When a HOBO problem is solved, the simulated bifurcation algorithm solves differential equations given by Equation (2), Equation (6), and Equation (7) using the symplectic Euler method.

When the symplectic Euler method is used, the differential equations given by Equation (2), Equation (3), and Equation (4) and the differential equations given by Equation (2), Equation (6), and Equation (7) can be written into a discrete recurrence relation as shown by Equation (8).

x_i(t+Δt)=x_i(t)+Dy_i(t)Δt
y_i(t+Δt)=y_i(t)+[{−D+p(t+Δt)−Kx_i²(t+Δt)}x_i(t+Δt)+f_i(t+Δt)]Δt
f_i(t+Δt)=−cz_i(t+Δt) (8)

t is time. At is a time step (unit time, time step size).

When the simulated bifurcation algorithm is executed, a digital computer or an electronic circuit such as an FPGA updates N variables x_iand N variables y_isequentially for each time step from the initial time and alternately between variables x_iand variables y_i, based on the algorithm in Equation (8). The digital computer or the electronic circuit such as an FPGA then binarizes the values of N variables x_iat the end time using the signum function and outputs the values of N spins.

Equation (8) is written using time t and time step Δt in order to show the correspondence with the differential equations. However, when the symplectic Euler method is performed by a digital computer or an electronic circuit such as an FPGA, the algorithm for computing Equation (8) may not necessarily include time t or time step Δt as explicit parameters. For example, when time step Δt is 1, the algorithm for computing Equation (8) need not include time step Δt. For example, when time t is not included as an explicit parameter, the algorithm for computing Equation (8) performs a process with x_i(t+Δt) as the updated value of x_i(t). More specifically, the algorithm for computing Equation (8) performs a process using “t” in Equation (8) as a parameter specifying a variable before update and “t+Δt” as a parameter specifying a variable after update. This is applicable to an improved algorithm of Equation (8) described below.

Improved Algorithm

The inventor has invented an algorithm (improved algorithm) by improving the simulated bifurcation algorithm. The improved algorithm further performs the following process when computing the differential equations given by Equation (8) using the symplectic Euler method.

More specifically, when the variable x_iis smaller than a predetermined first value after updating of the variable x_i, the improved algorithm changes the variable x_ito the first value and changes the variable y_ito a predetermined third value. When the variable x_iis greater than a second value greater than the first value after updating of the variable x_i, the improved algorithm changes the variable x_ito the second value and changes the variable y_ito the third value.

For example, the first value is −1, the second value is +1, and the third value is 0. In this case, when the variable x_iis smaller than −1 after updating of the variable x_i, the improved algorithm changes the variable x_ito −1 and changes the variable y_ito 0. When the variable x_iis greater than +1 after updating of the variable x_i, the improved algorithm changes the variable x_ito +1 and changes the variable y_ito 0.

Furthermore, after updating the variable x_iand the variable y_i, the improved algorithm adds an acceleration value calculated by a predetermined computation to the variable y_i. The acceleration value is a value, for example, obtained by multiplying the variable y_i(that is, the variable y_ibefore update) at the previous time by the product (rΔt) of a prescribed coefficient (r) and a unit time (Δt). The prescribed coefficient (r) is, for example, a predetermined constant. The prescribed coefficient may be a value determined for each updating based on a random number. When the first value is −1, the second value is +1, and the third value is 0, the prescribed coefficient is, for example, a value greater than 0 and smaller than 1.

“Kx²_i(t+Δt)” in Equation (8) is a term having the function of preventing divergence of the variable x_iduring calculation. Since the improved algorithm limits, for each updating, the variable x_iin a range from the first value (for example, −1) to the second value (for example, +1), both inclusive, the variable x_idoes not diverge during calculation. The improved algorithm therefore can compute a differential equation in which the term “Kx²_i(t+Δt)” in the differential equation given by Equation (8) is deleted.

Then, the improved algorithm is computed using Equation (9) instead of the algorithm in Equation (8).

x_i(t+Δt)=x_i(t)+Dy_i(t)Δt
y_i(t+Δt)=y_i(t)+[{−D+p(t+Δt)}x_i(t+Δt)+f_i(t+Δt)]Δt
f_i(t+Δt)=−cz_i(t+Δt)

In the case of a QUBO problem, z_i(t+Δt) in Equation (9) is given by Equation (10) below.

$\begin{matrix} z_{i} (t + Δ t) = - h_{i} α (t + Δ t) - \sum_{j = 1}^{N} J_{i, j} x_{j} (t + Δ t) & (10) \end{matrix}$

In the case of a HOBO problem, z_i(t+Δt) in Equation (9) is given by Equation (11) below.

$\begin{matrix} z_{i} (t + Δ t) = J_{i}^{(1)} α (t + Δ t) + \sum_{j = 1}^{N} J_{i, j}^{(2)} s_{j} (t + Δ t) + \sum_{j = 1}^{N} \sum_{k = 1}^{N} J_{i, j, k}^{(3)} s_{j} (t + Δ t) s_{k} (t + Δ t) + \dots & (11) \end{matrix}$

The improved algorithm in Equation (9) above is to solve the Hamiltonian equation in the same manner as Equation (8), and the variable y_icorresponds to the momentum. Therefore, the solution can be calculated stably using the symplectic Euler method, even without using a small value as the time step Δt.

The combinatorial optimization problem that can be solved by the improved algorithm is not limited to the Ising problem, and a common combinatorial optimization problem with binary variables can be solved. For example, the improved algorithm can be applied to a combinatorial optimization problem in which a variable of the objective function is a binary variable that takes one of a (first value) and b (second value) greater than a. When a solution of the objective function is calculated after a certain number of times of updating, a function f(x_i) whose range is binary, a or b, may be used instead of the signum function. The value of this function f(x_i) is determined based on the result of comparison of the value of the variable x_iwith a threshold value v (a<v<b). For example, if x_i<v, f(x_i)=a. If v<x_i, f(x_i)=b. For example, when x_i=v, f(x_i)=a or f(x_i)=b. The threshold value v may be, for example, (a+b)/2.

FIG. 1 is a diagram illustrating a bifurcation phenomenon of the variables x_iobtained when an optimization problem was solved by the simulated bifurcation algorithm. When an optimization problem is solved by the simulated bifurcation algorithm, a bifurcation phenomenon occurs, that is, a system with a single stable motion state changes to a system with two stable motion states as the parameters in the system change. As illustrated in FIG. 1, as the bifurcation phenomenon proceeds, the variables x_iare concentrated on the vicinity of −1 or +1 but also spread to a region smaller than −1 or a region greater than 1.

FIG. 2 is a diagram illustrating a bifurcation phenomenon of the variables x_iobtained when an optimization problem was solved by the improved algorithm. When an optimization problem is solved by the improved algorithm, a bifurcation phenomenon occurs, that is, a system with a single stable motion state changes to a system with two stable motion states as the parameters in the system change, in the same manner as the simulated bifurcation algorithm. However, the improved algorithm returns x_ismaller than −1 to −1 and returns x_igreater than +1 to +1, so that the variables x_ido not spread to a region smaller than −1 or a region greater than +1.

FIG. 3 and FIG. 4 are diagrams illustrating the result of the number of cuts with respect to the time step number when the maximum cut problem K2000 was solved 10000 times. FIG. 3 illustrates the mean value of the numbers of cuts. FIG. 4 illustrates the maximum value of the numbers of cuts. FIG. 5 illustrates the success probability that is the probability that the maximum cut value was obtained when the maximum cut problem K2000 was solved 10000 times using the improved algorithm.

FIG. 3, FIG. 4, and FIG. 5 illustrate the results of performing the improved algorithm where D is 1, K is 1, c is 0.9√(1999), r is 0.5, Δt is 1.1, and p is linearly increased from 0 to 1. The broken line Cmax in FIG. 3 and FIG. 4 indicates the maximum cut value of 33337 known in K2000. Referring to FIG. 3 and FIG. 4, it is understood that the improved algorithm comes close to the maximum cut value and can solve the optimization problem accurately. H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, K. Tatsumura, “High-performance combinatorial optimization based on classical mechanics”, Science Advances; 7, eabe7953”, February 2021 describes other methods (bSB, dSB) which are improved simulated bifurcation algorithms. FIG. 3, FIG. 4, and FIG. 5 also show the results of solving the maximum cut problem K2000 by bSB and dSB as comparative examples. The comparison with bSB and dSB has suggested that the improved algorithm provides equivalent or surpassing results. As illustrated in FIG. 5, it has also been found that the improved algorithm yields a higher success probability than bSB and dSB.

Functional Configuration

FIG. 6 is a diagram illustrating a functional configuration of a calculation device 10 according to the present embodiment.

The calculation device 10 solves an optimization problem using the improved algorithm. The calculation device 10 is implemented, for example, by an information processing device such as a computer, a computer system including a plurality of computers or servers mutually communicating through a network, or a PC cluster in which a plurality of computers work together to execute information processing. The calculation device 10 is implemented by electronic circuitry such as a CPU, a microprocessor, a GPU, an FPGA or an ASIC, or a combination thereof.

The calculation device 10 includes, as a functional configuration, an input module 12, an updating module 14, and an output module 16.

The input module 12 receives information (for example, N, J, h) for defining an objective function of a combinatorial optimization problem and information (x, Δt, T, p(t), α(t)) representing coefficients necessary for executing the improved algorithm from an external device and applies the received information to the updating module 14. The updating module 14 uses the improved algorithm and, for each of a plurality of elements each associated with the first variable (x_i) and the second variable (y_i), alternately updates the first variable (x_i) and the second variable (y_i) sequentially for each unit time (Δt) from an initial time (t=0) to an end time (t=T).

The output module 16 then outputs a solution to the combinatorial optimization problem based on the first variables (x_i) of the elements at the end time (t=T). For example, the output module 16 calculates, for each of the elements at the end time, a value of a discrete variable by binarizing the first variable (x_i) with a preset threshold value. The output module 16 then outputs the calculated values of discrete variables as a solution to the combinatorial optimization problem.

Process Flow

FIG. 7 is a diagram illustrating a first example of a process of the updating module 14. The updating module 14 performs a process, for example, according to the flow illustrated in FIG. 7.

First of all, at S101, the updating module 14 sets parameters. Specifically, the updating module 14 sets J, which is a matrix including N×N coupling coefficients, and h, which is an array including local magnetic field coefficients representing N local magnetic fields. When a HOBO problem is to be solved, the updating module 14 sets J⁽ⁿ⁾, which is an m^thrank tensor including Nm action coefficients, instead of J and h. In this case, (m) denotes the order of a variable of the objective function in the HOBO problem. The updating module 14 further sets the coefficient D, the coefficient c, Δt denoting the unit time, T denoting the end time, the function p(t), and the function α(t), where p(t) and α(t) are the increasing function that is 0 at t=initial time (for example, 0) and 1 at t=end time (T). The updating module 14 sets J and h in accordance with information received from the input module 12. The updating module 14 may set D, c, Δt, T, p(t), and α(t) in accordance with parameters received from the input module 12 or may set parameters determined in advance and unable to be changed.

Subsequently, at S102, the updating module 14 initializes the variables. Specifically, the updating module 14 initializes t that is a variable denoting time to the initial time (for example, 0). Furthermore, the updating module 14 substitutes an initial value received from the user, a predetermined fixed value, or a random number into each of N first variables (x₁(t) to x_N(t)) and N second variables (y₁(t) to y_N(t)).

Subsequently, the updating module 14 repeats the loop process between S103 and S118 until t becomes greater than T. In one loop process, the updating module 14 calculates N first variables (x₁(t+Δt) to x_N(t+Δt)) at the target time (t+Δt), based on N second variables (y_i(t+Δt) to y_N(t+Δt)) at the target time (t+Δt). In one loop process, the updating module 14 calculates N second variables (y_i(t+Δt) to y_N(t+Δt)) at the target time (t+Δt), based on N first variables (x₁(t) to x_N(t)) at the previous time (t) and N second variables (y₁(t) to y_N(t)) at the previous time (t).

The previous time (t) is the time a unit time (Δt) before the target time (t+Δt). That is, the updating module 14 repeats the loop process between S103 and S118 to update N first variables (x₁(t) to x_N(t)) and N second variables (y_i(t) to y_N(t)) sequentially for each unit time (Δt) from the initial time (t=0) to the end time (t=T).

At S104, the updating module 14 stores N second variables (y₁(t) to y_N(t)) at the previous time (t). The updating module 14 uses the stored N second variables (y_i(t) to y_N(t)) to calculate an acceleration value at S115 described later.

Subsequently, the updating module 14 repeats the loop process between S105 and S109 while incrementing i by one from i=1 to i=N, where i is an integer of 1 to N and an index representing the process target of N elements. Each of N elements is associated with the first variable (x_i(t)) and the second variable (y_i(t)). In the loop process between S105 and S109, the updating module 14 performs the process for the i^thelement of N elements as a target element.

At S106, the updating module 14 calculates an update value (z_i(t)), based on the first variable (x₁(t) to x_N(t)) at the previous time (t) for each of N elements and a predetermined action coefficient for each of sets of the target element and each of N elements. In the case of a QUBO problem, the action coefficient is the coupling coefficient included in J and the local magnetic field coefficient included in h. In the case of a HOBO problem, the action coefficient is included in J⁽ⁿ⁾.

In the case of a QUBO problem, the updating module 14 calculates Equation (12).

$\begin{matrix} z_{i} (t) = - h_{i} α (t) - \sum_{j = 1}^{N} J_{i, j} x_{j} (t) & (12) \end{matrix}$

In the case of a HOBO problem, the updating module 14 calculates Equation (13).

$\begin{matrix} z_{i} (t) = J_{i}^{(1)} α (t) + \sum_{j = 1}^{N} J_{i, j}^{(2)} x_{j} (t) + \sum_{j = 1}^{N} \sum_{k = 1}^{N} J_{i, j, k}^{(3)} x_{j} (t) x_{k} (t) + \dots & (13) \end{matrix}$

Subsequently, at S107, the updating module 14 calculates external force (f_i(t)) by multiplying the update value (z_i(t)) by the coefficient (c). Specifically, the updating module 14 calculates Equation (14).

f_i(t)=−cz_i(t) (14)

Subsequently, at S108, the updating module 14 calculates the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt) by adding a value obtained by multiplying a value based on the external force (f_i(t)) and the first variable (x_i(t)) of the target element at the previous time (t) by the unit time (Δt), to the second variable (y_i(t)) of the target element at the previous time (t). Specifically, the updating module 14 calculates Equation (15).

y_i(t+Δt)=y_i(t)+[{−D+p(t)}x(t)+f_i(t)]Δt (15)

For each of N elements, the updating module 14 updates the second variable (y_i(t+Δt)) at the target time (t+Δt), based on the first variable (x₁(t) to x_N(t)) at the previous time (t), by performing the loop process between S105 and S109 as described above N times. If the loop process between S105 and S109 is performed N times, the updating module 14 proceeds to S110.

Subsequently, the updating module 14 repeats the loop process between S110 and S116 while incrementing i by one from i=1 to i=N. In the loop process between S110 and S116, the updating module 14 performs the process for the i^thelement of N elements as a target element.

At S111, the updating module 14 calculates the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) by adding a value obtained by multiplying the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt) by the predetermined constant (D) and the unit time (Δt), to the first variable (x_i(t)) of the target element at the previous time (t). Specifically, the updating module 14 calculates Equation (16).

x_i(t+Δt)=x_i(t)+Dy_i(t+Δt)Δt (16)

Subsequently, at S112, the updating module 14 determines whether the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is in a range from a predetermined first value to a predetermined second value, both inclusive. In the present example, the first value is −1 and the second value is +1. For example, the second value is the unit amount of the first variable (x_i(t)) that is a continuous quantity. For example, the second value is a value with the reversed sign of the first value.

The updating module 14 proceeds to S115 (Yes at S112) if the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is in the range from the first value to the second value, both inclusive (Yes at S112). The updating module 14 proceeds to S113 if the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than the first value or greater than the second value (No at S112).

At S113, the updating module 14 performs a constraining process for the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt). Specifically, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than the first value, the updating module 14 changes the first variable (x_i(t+Δt)) to the first value. When the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is greater than the second value, the updating module 14 changes the first variable (x_i(t+Δt)) to the second value. In the present example, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than −1, the updating module 14 changes the first variable (x_i(t+Δt)) to −1. When the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is greater than +1, the updating module 14 changes the first variable (x_i(t+Δt)) to +1.

Subsequently to S113, at S114, the updating module 14 performs a constraining process for the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt). Specifically, the updating module 14 changes the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt) to a predetermined third value. The third value is, for example, a value greater than the first value and smaller than the second value. In the present example, the third value is 0. When S114 is finished, the updating module 14 proceeds to S115.

At S115, the updating module 14 adds an acceleration value calculated by a predetermined computation to the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt). For example, the acceleration value is a value obtained by multiplying the second variable (y_i(t)) at the previous time by the product (rΔt) of a prescribed coefficient (r) and a unit time (Δt). The prescribed coefficient (r) is, for example, a predetermined constant. For example, the prescribed coefficient (r) is a value greater than 0 and smaller than 1. In this case, the updating module 14 performs a computation given by Equation (17).

y_i(t+Δt)=y_i(t+Δt)+{rΔt×y_i(t)} (17)

The prescribed coefficient (r) may be any other value that does not depend on the time (t). For example, the prescribed coefficient (r) may be a value based on random numbers in a range greater than 0 and smaller than 1. When S115 is finished, the updating module 14 proceeds to S116.

The updating module 14 performs the following process by performing the loop process between S110 and S116 as described above N times. That is, for each of N elements, the updating module 14 updates the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt), based on the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt). Subsequently, for each of N elements, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than the first value, the updating module 14 changes the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) to the first value. Furthermore, for each of N elements, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is greater than the second value, the updating module 14 changes the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) to the second value. Furthermore, for each of N elements, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than the first value or greater than the second value, the updating module 14 changes the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt) to the third value. Furthermore, for each of N elements, the updating module 14 adds an acceleration value obtained by multiplying the second variable (y_i(t)) of the target element at the previous time (t) by the product (rΔt) of the prescribed coefficient (r) and the unit time (Δt), to the second variable (y_i(t+Δt)) of the target element at the target time (t+Δt).

If the loop process between S110 and S116 is performed N times, the updating module 14 proceeds to S117.

At S117, the updating module 14 updates the target time (t+Δt) by adding the unit time (Δt) to the previous time (t). At S118, the updating module 14 repeats the process from S104 to S117 until t exceeds the end time (T). When t becomes greater than the end time (T), the updating module 14 terminates this flow.

Then, for each of N elements, the updating module 14 calculates the value of the corresponding spin, in accordance with the sign of the first variable (x_i(T)) at the end time (t=T). For example, when the first variable (x_i(T)) at the end time (t=T) has a negative sign, the updating module 14 sets the corresponding spin to −1, and when positive, sets the corresponding spin to +1. Then, the updating module 14 outputs the calculated values of a plurality of spins or values obtained by converting the calculated values of a plurality of spins into discrete variables, as a solution to the combinatorial optimization problem.

By performing the process in accordance with the flowchart illustrated in FIG. 7, for each of N elements, the updating module 14 alternately updates the first variable (x_i(t+Δt)) and the second variable (y_i(t+Δt)), sequentially for each unit time from the initial time (t=0) to the end time (t=T). By performing the process in accordance with the flowchart in FIG. 7, in the updating process for each unit time, for each of N elements, the updating module 14 calculates the first variable (x_i(t+Δt)) at the target time (t+Δt) after calculating the second variable (y_i(t+Δt)) at the target time (t+Δt). By performing the process in accordance with the flowchart in FIG. 7, in the updating process for each unit time, for each of N elements, the updating module 14 calculates the second variable (y_i(t+Δt)) at the target time (t+Δt) before updating the first variable (x_i(t+Δt)) at the target time (t+Δt) and adds an acceleration value to the second variable (y_i(t+Δt)) at the target time (t+Δt) after updating the first variable (x_i(t+Δt)) at the target time (t+Δt).

As described above, the updating module 14 can calculate N first variables (x₁(t) to x_N(t)) and N second variables (y₁(t) to y_N(t)) at the end time (t=T) by performing a computation in accordance with the improved algorithm using the symplectic Euler method.

FIG. 8 is a diagram illustrating a second example of a process of the updating module 14. When an optimization problem is solved using the improved algorithm, the updating module 14 may perform a process according to the flow illustrated in FIG. 8, instead of the flow illustrated in FIG. 7. In the flowchart in FIG. 8, the step that performs substantially the same process as the process in the flowchart illustrated in FIG. 7 is denoted by the same step number.

First of all, at S101 and S102, the updating module 14 performs a process similar to that in the first example in FIG. 7. Subsequently, the updating module 14 performs the loop process between S103 and S118, in the same manner as in the first example in FIG. 7.

At S104, the updating module 14 performs a process similar to that in the first example illustrated in FIG. 7.

At S111, the updating module 14 calculates the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) by adding a value obtained by multiplying the second variable (y_i(t)) of the target element at the previous time (t) by the predetermined constant (D) and the unit time (Δt), to the first variable (x_i(t)) of the target element at the previous time (t). Specifically, the updating module 14 calculates Equation (18).

x_i(t+Δt)=x_i(t)+Dy_i(t)Δt (18)

Subsequently, at S112, the updating module 14 determines whether the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is in a range from a predetermined first value to a predetermined second value, both inclusive. The updating module 14 proceeds to S115 (Yes at S112) if the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is in the range from the first value to the second value, both inclusive (Yes at S112). The updating module 14 proceeds to S113 if the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than the first value or greater than the second value (No at S112).

At S113, the updating module 14 performs a process similar to that in the first example illustrated in FIG. 7. Subsequently to S113, at S114, the updating module 14 performs the constraining process for the second variable (y_i(t)) of the target element at the previous time (t). Specifically, the updating module 14 changes the second variable (y_i(t)) of the target element at the previous time (t) to the third value. When S114 is finished, the updating module 14 proceeds to S115.

At S115, the updating module 14 adds an acceleration value calculated by a predetermined computation to the second variable (y_i(t)) of the target element at the previous time (t). For example, the acceleration value is a value obtained by multiplying the second variable (y′_i(t)) at the previous time before performing the constraining process by the product (rΔt) of the prescribed coefficient (r) and the unit time (Δt). The updating module 14 stores the second variable (y′_i(t)) at the previous time before performing the constraining process in advance at S104. For example, the updating module 14 performs a computation given by Equation (19).

y_i(t)=y_i(t)+{rΔt×y′_i(t)} (19)

When S115 is finished, the updating module 14 proceeds to S116. The updating module 14 performs the following process by performing the loop process between S110 and S116 as described above N times.

That is, for each of N elements, the updating module 14 updates the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt), based on the second variable (y_i(t)) of the target element at the previous time (t). Subsequently, for each of N elements, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than the first value, the updating module 14 changes the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) to the first value. Furthermore, for each of N elements, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is greater than the second value, the updating module 14 changes the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) to the second value. Furthermore, for each of N elements, when the first variable (x_i(t+Δt)) of the target element at the target time (t+Δt) is smaller than the first value or greater than the second value, the updating module 14 changes the second variable (y_i(t)) of the target element at the previous time (t) to the third value. Furthermore, for each of N elements, the updating module 14 adds an acceleration value obtained by multiplying the second variable (y_i(t)) of the target element at the previous time (t) before performing the constraining process by the product (rΔt) of the prescribed coefficient (r) and the unit time (Δt), to the second variable (y_i(t)) of the target element at the previous time (t).

If the loop process between S110 and S116 is performed N times, the updating module 14 proceeds to S105.

Subsequently, the updating module 14 repeats the loop process between S105 and S109 while incrementing i by one from i=1 to i=N.

At S106, the updating module 14 calculates an update value (z_i(t+Δt)), based on the first variable (x₁(t+Δt) to x_N(t+Δt)) at the target time (t+Δt) for each of N elements and a predetermined action coefficient for each of sets of the target element and each of N elements.

In the case of a QUBO problem, the updating module 14 calculates Equation (20).

$\begin{matrix} z_{i} (t + Δ t) = - h_{i} α (t + Δ t) - \sum_{j = 1}^{N} J_{i, j} x_{j} (t + Δ t) & (20) \end{matrix}$

In the case of a HOBO problem, the updating module 14 calculates Equation (21).

$\begin{matrix} z_{i} (t + Δ t) = J_{i}^{(1)} α (t + Δ t) + \sum_{j = 1}^{N} J_{i, j}^{(2)} x_{j} (t + Δ t) + \sum_{j = 1}^{N} \sum_{k = 1}^{N} J_{i, j, k}^{(3)} x_{j} (t + Δ t) x_{k} (t + Δ t) + \dots & (21) \end{matrix}$

Subsequently, at S107, the updating module 14 calculates the external force (f_i(t+Δt)) by multiplying the update value (z_i(t+Δt)) by the coefficient (c). Specifically, the updating module 14 calculates Equation (22).

f_i(t+Δt)=−cz_i(t+Δt) (22)

Subsequently, at S108, the updating module 14 calculates the second variable (y_i(t+A)) of the target element at the target time (t+Δt) by adding a value obtained by multiplying a value based on the external force (f_i(t)) and the first variable (x_i(t+A)) of the target element at the target time (t+A) by the unit time (Δt), to the second variable (y_i(t)) of the target element at the previous time (t). Specifically, the updating module 14 calculates Equation (23).

y_i(t+Δt)=y_i(t)+[{−D+p(t+Δt)}x_i(t+Δt)+f_i(t+Δt)]Δt (23)

For each of N elements, the updating module 14 updates the second variable (y_i(t+Δt)) at the target time (t+Δt), based on the first variable (x₁(t+Δt) to x_N(t+Δt)) at the target time (t+Δt), by performing the loop process between S105 and S109 as described above N times.

If the loop process between S105 and S109 is performed N times, the updating module 14 proceeds to S117.

Then, for each of N elements, the updating module 14 calculates the value of the corresponding spin, in accordance with the sign of the first variable (x_i(T)) at the end time (t=T).

By performing the process in accordance with the flowchart illustrated in FIG. 8, for each of N elements, the updating module 14 alternately updates the first variable (x_i(t+Δt)) and the second variable (y_i(t+Δt)), sequentially for each unit time from the initial time (t=0) to the end time (t=T). By performing the process in accordance with the flowchart in FIG. 8, in the updating process for each unit time, for each of N elements, the updating module 14 calculates the second variable (y_i(t+Δt)) at the target time (t+Δt) after calculating the first variable (x_i(t+Δt)) at the target time (t+Δt). By performing the process in accordance with the flowchart in FIG. 8, in the updating process for each unit time, for each of N elements, the updating module 14 adds the acceleration value to the second variable (y_i(t)) at the previous time (t) after updating the first variable (x_i(t+Δt)) and before updating the second variable (y_i(t+Δt)) at the target time (t+Δt).

The updating module 14 performs the process at S110 to S116 for each of N elements. The updating module 14 may allow M₁(where M₁is an integer of 2 through N, both inclusive) computation processing devices such as processing circuits, processors, computers, or servers to perform these processes in parallel. For example, M₁is a divisor of N. In this case, each of N elements corresponds to any one of M₁computation processing devices. Then, each of M₁computation processing devices performs the process at S110 to S116 for the corresponding element among N elements. By performing such a process, the updating module 14 allows a plurality of circuits such as FPGAs to perform these processes in a distributed manner and thereby obtains a processing result at high speed.

Furthermore, the updating module 14 also performs the process at S105 to S109 for each of N elements. The updating module 14 may allow M₂computation processing devices such as processing circuits, processors, computers, or servers to perform these processes in parallel. For example, M₂is a divisor of N. In this case, each of N elements corresponds to any one of M₂computation processing devices. Then, each of M₂computation processing devices performs the process at S105 to S109 for the corresponding element among N elements. By performing such a process, the updating module 14 allows a plurality of circuits such as FPGAs to perform these processes in a distributed manner and thereby obtains a processing result at high speed.

Computing Circuit

FIG. 9 is a diagram illustrating a block configuration of the updating module 14 implemented in circuitry. The updating module 14 implemented in circuitry is configured as illustrated in FIG. 9 as an example. In the description of FIG. 9, the previous time is denoted by t₁. The target time is denoted by t₂.

The updating module 14 illustrated in FIG. 9 executes the improved algorithm given by Equation (9) and Equation (10). The updating module 14 includes an X memory 66, a Y memory 67, an action computation circuit 68, an updating circuit 69, and a control circuit 70.

The X memory 66 stores N first variables (x_i(t₁)) at the previous time (t₁). The N first variables (x_i(t₁)) at the previous time (t₁) stored in the X memory 66 are overwritten with updating of the time. That is, when N first variables (x_i(t₂)) at the target time (t₂) are calculated, the calculated N first variables (x_i(t₂)) at the target time (t₂) are written into the X memory 66 as new N first variables (x_i(t₁)) at the previous time (t₁). Prior to computation, the control circuit 70 writes N first variables x_iat the initial time into the X memory 66.

The Y memory 67 stores N second variables (y_i(t₁)) at the previous time (t₁). N second variables (y_i(t₁)) at the previous time (t₁) stored in the Y memory 67 are overwritten with updating of the time. That is, when N second variables (y_i(t₂)) at the target time (t₂) are calculated, the calculated N second variables (y_i(t₂)) at the target time (t₂) are written into the Y memory 67 as new N second variables (y_i(t₁)) at the previous time (t₁). Prior to computation, the control circuit 70 writes N second variables y_iat the initial time into the Y memory 67.

The action computation circuit 68 acquires N first variables (x_j(t₁)) at the previous time (t₁) from the X memory 66. Then, for each of N elements, the action computation circuit 68 calculates the update value (z_i(t₁)).

For each of N elements, the updating circuit 69 acquires the update value (z_i(t₁)) from the action computation circuit 68. Furthermore, for each of N elements, the updating circuit 69 acquires the first variable (x_i(t₁)) at the previous time (t₁) from the X memory 66 and acquires the second variable (y_i(t₁)) at the previous time (t₁) from the Y memory 67. Then, for each of N elements, the updating circuit 69 calculates the first variable (x_i(t₂)) at the target time (t₂) and overwrites the first variable (x_i(t₁)) at the previous time (t₁) stored in the X memory 66. In addition, for each of N elements, the updating circuit 69 calculates the second variable (y_i(t₂)) at the target time (t₂) and overwrites the second variable (y_i(t₁)) at the previous time (t₁) stored in the Y memory 67.

The control circuit 70 sequentially updates the target time (t₂) for each unit time (Δt) to allow the action computation circuit 68 and the updating circuit 69 to sequentially calculate the first variable (x_i(t)) and the second variable (y_i(t)) for each unit time (Δt).

Furthermore, the control circuit 70 generates an index (i) by incrementing from 1 to N and allows the action computation circuit 68 and the updating circuit 69 to calculate the first variable (x_i(t₂)) at the target time (t₂) and the second variable (y_i(t₂)) at the target time (t₂) corresponding to each of N elements, in the order of index. The action computation circuit 68 and the updating circuit 69 may calculate a plurality of first variables (x_i(t₂)) and a plurality of second variables (y_i(t₂)) corresponding to a plurality of indices, in parallel.

The action computation circuit 68 includes a J memory 71, an H memory 72, a matrix computation circuit 73, an a function circuit 74, and a first adder circuit 75.

The J memory 71 stores an N×N matrix including (N×N) coupling coefficients. J_i,jdenotes the coupling coefficient at the i^throw and the j^thcolumn included in the matrix. J_i,jdenotes the coupling coefficient between the i^thspin and the j^thspin in the Ising model corresponding to a combinatorial optimization problem. Prior to computation, the control circuit 70 writes a matrix generated by a user or the like in advance into the J memory 71.

The H memory 72 stores an array including N local magnetic field coefficients. h_idenotes the i^thlocal magnetic field coefficient included in the array. h_idenotes the local magnetic field acting on the i^thspin in the Ising model corresponding to a combinatorial optimization problem. Prior to computation, the control circuit 70 writes an array generated by a user or the like in advance into the H memory 72.

The matrix computation circuit 73 acquires N first variables (x_j(t₁)) at the previous time (t₁) from the X memory 66. The matrix computation circuit 73 acquires, for each of N elements, N coupling coefficients J_i,jincluded in a target row from the J memory 71. Then, the matrix computation circuit 73 performs, for each of N elements, a product-sum operation of N first variables (x_j(t₁)) at the previous time (t₁) and N coupling coefficients J_i,jincluded in the target row.

The α function circuit 74 acquires, for each of N elements, the target local magnetic field coefficient h_ifrom the H memory 72. The α function circuit 74 performs, for each of N elements, a computation {−h_iα(t₁)}, where α(t) is a preset function.

The first adder circuit 75 adds, for each of N elements, the result of the product-sum operation by the matrix computation circuit 73 to the computation result by the α function circuit 74. With this computation, the first adder circuit 75 outputs, for each of N elements, the update value (z_i(t₁)) at the previous time (t₁) given by Equation (24).

$\begin{matrix} z_{i} (t_{1}) = - h_{i} α (t_{1}) - \sum_{j = 1}^{N} J_{i, j} x_{j} (t_{1}) & (24) \end{matrix}$

The updating circuit 69 includes a first multiplier circuit 79, a P function circuit 80, a second multiplier circuit 81, a second adder circuit 82, a third multiplier circuit 83, a third adder circuit 84, a before-constraint Y memory 85, a fourth multiplier circuit 86, a fourth adder circuit 87, a determination circuit 91, an X constraint circuit 92, a Y constraint circuit 93, an acceleration value calculation circuit 94, and an acceleration value addition circuit 95.

The first multiplier circuit 79 multiplies, for each of N elements, the update value (z_i(t₁)) by a coefficient −c. The P function circuit 80 performs, for each of N elements, a computation of {−D+p(t₁)}. The second multiplier circuit 81 acquires, for each of N elements, the first variable (x_i(t₁)) at the previous time (t₁) from the X memory 66. Then, the second multiplier circuit 81 multiplies, for each of N elements, the first variable (x_i(t₁)) at the previous time (t₁) by the computation result in the P function circuit 80.

The second adder circuit 82 adds, for each of N elements, the computation result in the first multiplier circuit 79 to the computation result in the second multiplier circuit 81. The third multiplier circuit 83 multiplies the computation result in the second adder circuit 82 by the unit time Δt.

The third adder circuit 84 acquires, for each of N elements, the second variable (y_i(t₁)) at the previous time (t₁) from the Y memory 67. The third adder circuit 84 adds, for each of N elements, the second variable (y_i(t₁)) at the previous time (t₁) to the computation result in the third multiplier circuit 83. With this computation, the third adder circuit 84 outputs, for each of N elements, the second variable (y_i(t₂)) at the target time (t₂) given by Equation (25).

y_i(t₂)=y_i(t₁)+[{−D+p(t₁)}x_i(t)−cz_i(t₁)]Δt (25)

Then, the third adder circuit 84 writes the calculated second variable (y_i(t₂)) at the target time (t₂) for each of N elements into the before-constraint Y memory 85. The before-constraint Y memory 85 stores N second variables (y_i(t₂)) at the target time (t₂) before constraint by the Y constraint circuit 93.

The fourth multiplier circuit 86 acquires, for each of N elements, the second variable (y_i(t₂)) before constraint at the target time (t₂) from the before-constraint Y memory 85. The fourth multiplier circuit 86 multiplies, for each of N elements, the second variable (y_i(t₂)) before constraint at the target time (t₂) by {DΔt}.

The fourth adder circuit 87 acquires, for each of N elements, the first variable (x_i(t₁)) at the previous time (t₁) from the X memory 66. The fourth adder circuit 87 adds, for each of N elements, the first variable (x_i(t₁)) at the previous time (t₁) to the computation result in the fourth multiplier circuit 86. With this computation, the fourth adder circuit 87 outputs, for each of N elements, the first variable (x_i(t₂)) at the target time (t₂) given by Equation (26).

x_i(t₂)=x_i(t₁)+Dy_i(t₂)Δt (26)

The determination circuit 91 determines, for each of N elements, whether the first variable (x_i(t₂)) at the target time (t₂) calculated by the fourth adder circuit 87 is in a range from a predetermined first value to a predetermined second value, both inclusive. For example, the first value is −1. The second value is +1. When the first variable (x_i(t₂)) at the target time (t₂) is smaller than the first value or greater than the second value, the determination circuit 91 applies an enable signal (EN) to the X constraint circuit 92 and the Y constraint circuit 93.

The X constraint circuit 92 receives, for each of N elements, the first variable (x_i(t₂)) at the target time (t₂) calculated by the fourth adder circuit 87. For each of N elements, when an enable signal (EN) is not received from the determination circuit 91, the X constraint circuit 92 writes the first variable (x_i(t₂)) at the target time (t₂) calculated by the fourth adder circuit 87 as it is into the X memory 66.

For each of N elements, when an enable signal (EN) is received from the determination circuit 91, the X constraint circuit 92 performs the constraining process for the first variable (x_i(t₂)) at the target time (t₂) calculated by the fourth adder circuit 87 and writes the first variable (x_i(t₂)) subjected to the constraining process into the X memory 66.

When the first variable (x_i(t₂)) at the target time (t₂) is smaller than the first value (for example, −1), the X constraint circuit 92 changes the first variable (x_i(t₂)) at the target time (t₂) to the first value (for example, −1), as the constraining process. When the first variable (x_i(t₂)) at the target time (t₂) is greater than the second value (for example, +1), the X constraint circuit 92 changes the first variable (x_i(t₂)) at the target time (t₂) to the second value (for example, +1), as the constraining process.

The Y constraint circuit 93 acquires, for each of N elements, the second variable (y_i(t₂)) at the target time (t₂) from the before-constraint Y memory 85. For each of N elements, when an enable signal (EN) is not received from the determination circuit 91, the Y constraint circuit 93 outputs the second variable (y_i(t₂)) at the target time (t₂) acquired from the before-constraint Y memory 85 as it is.

For each of N elements, when an enable signal (EN) is received from the determination circuit 91, the Y constraint circuit 93 performs the constraining process for the second variable (y_i(t₂)) at the target time (t₂) acquired from the before-constraint Y memory 85 and outputs the second variable (y_i(t₂)) subjected to the constraining process.

Here, the Y constraint circuit 93 changes the second variable (y_i(t₂)) at the target time (t₂) to the third value (for example, 0), as the constraining process. The Y constraint circuit 93 may change the second variable (y_i(t₂)) to a random number from the first value to the second value, both inclusive, as the constraining process.

The acceleration value calculation circuit 94 acquires, for each of N elements, the second variable (y_i(t₁)) at the previous time (t₁) from the Y memory 67. Furthermore, the acceleration value calculation circuit 94 receives the prescribed coefficient (r) from the control circuit 70 prior to computation. For example, the prescribed coefficient (r) is a value greater than 0 and smaller than 1. The acceleration value calculation circuit 94 then calculates an acceleration value, for each of N elements, by multiplying the second variable (y_i(t₁)) at the previous time (t₁) by the product (rΔt) of the prescribed coefficient (r) and the unit time (Δt).

The acceleration value addition circuit 95 acquires, for each of N elements, the second variable (y_i(t₂)) at the target time (t₂) from the Y constraint circuit 93. Furthermore, the acceleration value addition circuit 95 acquires, for each of N elements, the acceleration value (rΔt×(y_i(t₁))) from the acceleration value calculation circuit 94. The acceleration value addition circuit 95 adds, for each of N elements, the acceleration value (rΔt×(y_i(t₁))) to the second variable (y_it₂)) at the target time (t₂). The acceleration value addition circuit 95 then writes the second variable (y_i(t₂)) at the target time (t₂) added with the acceleration value (rΔt×(y_i(t₁))) into the Y memory 67.

As described above, the updating module 14 can calculate N first variables (x_i(T)) and N second variables (y_i(T)) at the end time (T) by executing the improved algorithm in Equation (9).

The updating module 14 described above calculates, for each unit time, the first variable (x_i(t₂)) at the target time (t₂) after calculating the second variable (y_i(t₂)) at the target time (t₂). Instead of this, the updating module 14 may calculate, for each unit time, the second variable (y_i(t₂)) at the target time (t₂) after calculating the first variable (x_i(t₂)) at the target time (t₂).

In this case, the before-constraint Y memory 85 stores N second variables (y_i(t₁)) at the previous time (t₁) before constraint by the Y constraint circuit 93. Then, the fourth adder circuit 87 outputs, for each of N elements, the first variable (x_i(t₂)) at the target time (t₂) given by Equation (27).

x_i(t₂)=x_i(t₁)+Dy_i(t₁)Δt (27)

The matrix computation circuit 73 acquires N first variables (x_j(t₂)) at the target time (t₂) from the X memory 66. The first adder circuit 75 then outputs, for each of N elements, the update value (z_i(t₂)) at the target time (t₂) given by Equation (28).

$\begin{matrix} z_{i} (t_{2}) = - h_{i} α (t_{2}) - \sum_{j = 1}^{N} J_{i, j} x_{j} (t_{2}) & (28) \end{matrix}$

The first multiplier circuit 79 multiplies, for each of N elements, the update value (z_i(t₂)) at the target time (t₂) by −c. The P function circuit 80 performs, for each of N elements, a computation of {−D+p(t₂)}. The second multiplier circuit 81 multiplies, for each of N elements, the first variable (x_i(t₂)) at the previous time (t₂) by the computation result in the P function circuit 80.

The third adder circuit 84 then outputs, for each of N elements, the second variable (y_i(t₂)) at the target time (t₂) given by Equation (29).

y_i(t₂)=y_i(t₁)+[{−D+p(t₂)}x_i(t₂)−cz_i(t₂)]Δt (29)

Even with such a process, the updating module 14 can calculate N first variables (x_i(T)) and N second variables (y_i(T)) at the end time (T) by executing the improved algorithm.

The updating circuit 69 performs, for each of N elements, an X updating process of updating the first variable (x_i(t₁)) based on the second variable (y_i(t₁)), an X constraining process of changing the first variable (x_i(t₂)) to the first value or the second value, a Y constraining process of changing the second variable (y_i(t₂)) to the third value, and an acceleration value adding process of adding the acceleration value (rΔt×(y_i(t₁))) to the second variable (y_i(t₂)). The updating circuit 69 may perform these processes with M₁processing circuits in parallel. In this case, each of N elements corresponds to any one of M₁processing circuits. Then, each of M₁processing circuits performs, for the corresponding element among N elements, the X updating process, the X constraining process, the Y constraining process, and the acceleration value adding process. By performing such a process, the updating circuit 69 allows a plurality of processing circuits to perform these processes in a distributed manner and thereby obtains a processing result at high speed.

Furthermore, the action computation circuit 68 calculates, for each of N elements, the update value (z_i(t₁)). The action computation circuit 68 may allow M₂processing circuits to perform the process of calculating the update value (z_i(t₁)) in parallel. In this case, each of N elements corresponds to any one of M₂processing circuits. Then, each of M₂processing circuits performs the process of calculating the update value (z_i(t₁)) for the corresponding element among N elements. When the action computation circuit 68 is parallelized, the degree of parallelism of the matrix computation circuit 73, the degree of parallelism of the α function circuit 74, and the degree of parallelism of the first adder circuit 75 may vary, or only the matrix computation circuit 73 may be parallelized. By performing such a process, the action computation circuit 68 allows a plurality of processing circuits to perform the process of calculating the update value (z_i(t₁)) in a distributed manner and thereby obtains a processing result at high speed. In particular, with the parallelization of the matrix computation circuit 73, the action computation circuit 68 can efficiently perform an enormous amount of multiplication and addition processes involved with the matrix computation.

System Configuration

FIG. 10 is a diagram illustrating a configuration of an information processing system 100. The improved algorithm can be executed by, for example, the information processing system 100 as illustrated in FIG. 10. The information processing system 100 thus can solve a large-scale combinatorial optimization problem 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 devices), a plurality of cables 104 (104a to 104c), and a switch 105. A terminal device 106 can communicate with the information processing system 100. The management server 101, the calculation servers 103 (103a to 103c), and the terminal device 106 perform data communication with each other through the network 102. The network 102 is, for example, the Internet in which a plurality of computer networks are connected to each other. The network 102 may be a communication medium either wired or wireless or a combination thereof.

The calculation servers 103 (103a to 103c) are connected to the switch 105 through the cables 104 (104a to 104c). The cables 104 (104a to 104c) and the switch 105 form an interconnection between the calculation servers 103. The calculation servers 103 (103a to 103c) can perform data communication with each other through the interconnection. The switch 105 is, for example, an InfiniBand switch. The cables 104a to 104c are, for example, InfiniBand cables. However, the switch 105 and the cables 104 may be wired LAN switch/cables. Any communication standards and communication protocol may be used in the cables 104 and the switch 105. The terminal device 106 is, for example, a notebook PC, a desktop PC, a smartphone, a tablet, or a vehicle-mounted terminal device.

The process of solving a combinatorial optimization problem using the improved algorithm can be performed in parallel and in a distributed manner. Therefore, the calculation servers 103 (103a to 103c) and/or the processors of the calculation servers 103 (103a to 103c) may share and execute some of the steps of some calculation processes in the process of solving a combinatorial optimization problem using the improved algorithm or may perform the same calculation process for different variables in parallel. In this case, the management server 101, for example, converts a combinatorial optimization problem input by a user into a format that can be processed by the calculation servers 103 and allows the calculation servers 103 to perform a process. The management server 101 then acquires the calculation results from the calculation servers 103 and converts the consolidated calculation result into the solution to the combinatorial optimization problem.

FIG. 10 illustrates three calculation servers 103 (103a to 103c). However, the number of calculation servers 103 included in the information processing system 100 is not limited to three. For example, the information processing system 100 may include one, two, four, or more calculation servers 103. The information processing system 100 may perform the process of solving a combinatorial optimization problem using some of the included calculation servers 103. The calculation server 103 may be any kind of information processing device. For example, the calculation server 103 may be a server installed in a datacenter or may be a desktop PC installed in an office. The calculation server 103 may be computers of different kinds installed at different locations. For example, the calculation server 103 may be a general-purpose computer, or may be a dedicated electronic circuit or a combination thereof.

FIG. 11 is a diagram illustrating a configuration of the management server 101. The management server 101 is, for example, a computer including a central processing unit (CPU) and a 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. The processor 110, the storage unit 114, the communication circuit 115, the input circuit 116, and the output circuit 117 are connected to each other through a bus 120. The processor 110 includes, as an internal functional configuration, a management module 111, a conversion module 112, and a control module 113.

The processor 110 is an electronic circuit that performs a computation and controls the management server 101. The processor 110 may be, for example, a CPU, a microprocessor, an ASIC, an FPGA, a PLD, or a combination thereof. The management module 111 provides an interface for operating the management server 101 through the user's terminal device 106. The interface provided by the management module 111 is, for example, an API, a CLI, or a webpage. For example, the user inputs information on a combinatorial optimization problem or views and/or downloads the calculated solution to the combinatorial optimization problem through the management module 111. The conversion module 112 receives parameters for a combinatorial optimization problem and converts the input parameters into a form that can be processed by each calculation server 103. The control module 113 transmits a control command to each calculation server 103. After the control module 113 acquires a calculation result from each calculation server 103, the conversion module 112 consolidates a plurality of calculation results, converts the consolidated calculation result into a solution to the combinatorial optimization problem, and outputs the solution to the combinatorial optimization problem.

The storage unit 114 stores a computer program for the management server 101 and a variety of data including data necessary for executing the computer program and data generated by the computer program. The computer program includes both an OS and an application. The storage unit 114 may be a volatile memory, a nonvolatile memory, or a combination thereof. The volatile memory is, for example, a DRAM or an SRAM. The nonvolatile memory is, for example, a NAND flash memory, an NOR flash memory, an ReRAM, or an MRAM. Alternatively, the storage unit 114 may be a hard disk, an optical disk, a magnetic tape, or an external storage device.

The communication circuit 115 transmits/receives data to/from the devices connected to the network 102. The communication circuit 115 is, for example, a network interface card (NIC) for a wired LAN. However, the communication circuit 115 may be any other kind 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 or PCI-Express as an external port. An operation 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, but not limited to, a keyboard, a mouse, a touch panel, or a voice recognition device. The output circuit 117 implements data output from the management server 101. The output circuit 117 includes, for example, HDMI (registered trademark) or DisplayPort as an external port. For example, a display device 119 is connected to the output circuit 117. The display device 119 is, for example, but not limited to, a liquid crystal display (LCD), an organic electroluminescent (EL) display, or a projector.

An administrator can perform maintenance of the management server 101, using the operation device 118 and the display device 119. The operation device 118 and the display device 119 may be built in the management server 101. The operation device 118 and the display device 119 are not necessarily connected to the management server 101. For example, the administrator may perform maintenance of the management server 101, using a terminal device capable of communicating with the network 102.

FIG. 12 is a diagram illustrating data stored in the storage unit 114 of the management server 101. The storage unit 114 stores, for example, problem data 114A, calculation data 114B, a management program 114C, a conversion program 114D, and a control program 114E. For example, the problem data 114A includes data of a combinatorial optimization problem. For example, the calculation data 114B includes a calculation result collected from each calculation server 103. For example, the management program 114C is a computer program that implements the function of the management module 111. For example, the conversion program 114D is a computer program that implements the function of the conversion module 112. For example, the control program 114E is a computer program that implements the function of the control module 113.

FIG. 13 is a diagram illustrating a configuration of the calculation server 103a. Other calculation servers 103 may have a configuration similar to that of the calculation server 103a or may have a configuration different from that of 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 adaptor 135. The communication circuit 131, the shared memory 132, the processors 133a to 133d, the storage 134, and the host bus adaptor 135 are connected to each other through a bus 136.

The communication circuit 131 transmits/receives data to/from the devices connected to the network 102. The communication circuit 131 is, for example, a network interface card (NIC) for a wired LAN. However, the communication circuit 131 may be any other kind of communication circuit, such as a wireless LAN. The shared memory 132 is a memory accessible by the processors 133a to 133d. The shared memory 132 is, for example, a volatile memory such as a DRAM and an SRAM. However, the shared memory 132 may include any other kind of memory such as a nonvolatile memory. The processors 133a to 133d share data through the shared memory 132. The shared memory 132 is not necessarily configured with all of the memories of the calculation server 103a. For example, some memories of the calculation servers 103a may be local memories accessed only by any one of the processors 133a to 133d.

The processors 133a to 133d are electronic circuits that perform a calculation process. The processors 133a to 133d may be, for example, CPUs, GPUs, FPGAs, or ASICs, or may be a combination thereof. The processors 133a to 133d may be CPU cores or CPU threads. When the processors 133a to 133d are CPUs, the number of sockets included in the calculation server 103a is not limited. The processors 133a to 133d may be connected to any other components of the calculation server 103a through a bus such as PCI Express.

In the example in FIG. 13, the calculation server 103a includes four processors 133a to 133d. However, the number of processors included in one calculation server 103a is not limited to four.

The storage 134 stores a computer program for the calculation server 103a and a variety of data including data necessary for executing the computer program and data generated by the computer program. As used herein the computer program includes both an OS and an application. The storage 134 may be a volatile memory and a nonvolatile memory, or a combination thereof. The volatile memory is, for example, a DRAM or an SRAM. The nonvolatile memory is, for example, a NAND flash memory, an NOR flash memory, an ReRAM, or an MRAM. The storage 134 may include a hard disk, an optical disk, a magnetic tape, or an external storage device.

The host bus adaptor 135 implements data communication with other calculation servers 103. The host bus adaptor 135 is connected to the switch 105 through the cable 104a. The host bus adaptor 135 is, for example, a host channel adaptor (HCA). The host bus adaptor 135, the cable 104a, and the switch 105 form an interconnection that can achieve a high throughput and thereby can improve the speed of parallel calculation processing.

FIG. 14 is a diagram illustrating data stored in the storage 134. The storage 134 stores, for example, calculation data 134A, a calculation program 134B, and a control program 134C. The calculation data 134A includes data in progress of calculation and a calculation result of the calculation server 103a. At least part of the calculation data 134A may be stored in a different storage level, such as the shared memory 132, a cache of a processor, or a register of a processor. The calculation program 134B is a computer program that implements a calculation process in each processor 133 and a storage process of data into the shared memory 132 and the storage 134. The control program 134C is a computer program that controls the calculation server 103a based on a command transmitted from the control module 113 of the management server 101 and transmits the calculation result of the calculation server 103a to the management server 101.

Such a calculation server 103a executes a computer program for the processors 133a to 133d to solve a combinatorial optimization problem. This computer program causes the calculation server 103a to function as the input module 12, the updating module 14, and the output module 16.

The information processing system 100 having such a configuration can be used as a PC cluster. The PC cluster is a system in which a plurality of computers are connected to achieve calculation performance unachievable by one computer.

The information processing system 100 can perform parallel calculation even in a configuration in which memories are distributed over a plurality of calculation servers 103, for example, using a message passing interface (MPI).

The information processing system 100 may be a plurality of GPUs connected through a high-speed link. In this case, each of the GPUs performs a process similar to that of the calculation server 103.

The present invention is not limited to the foregoing embodiments as they are and, in the implementation phase, can be embodied with constituent elements modified without departing from the spirit of the invention. A plurality of constituent elements disclosed in the foregoing embodiments may be combined as appropriate to form a variety of inventions. For example, some of all of the constituent elements described in the embodiments may be eliminated. The constituent elements in different embodiments may be combined as appropriate.

While certain embodiments have been described, these embodiments have been presented by way of example only, and are not intended to limit the scope of the inventions. Indeed, the novel embodiments described herein may be embodied in a variety of other forms; furthermore, various omissions, substitutions and changes in the form of the embodiments described herein may be made without departing from the spirit of the inventions. The accompanying claims and their equivalents are intended to cover such forms or modifications as would fall within the scope and spirit of the inventions.

TECHNICAL IDEAS

The foregoing embodiments can be summarized in the following technical ideas.

Technical Idea 1

A calculation device configured to solve a combinatorial optimization problem, the calculation device including:

a memory; and

one or more processors coupled to the memory and configured to:

- alternately update, for a plurality of elements each associated with a first variable and a second variable, the first variable and the second variable, sequentially for unit times from an initial time to an end time; and
- output a solution to the combinatorial optimization problem based on the first variables of the plurality of elements at the end time, wherein

the plurality of elements correspond to a plurality of discrete variables of the combinatorial optimization problem,

the first variables and the second variables are each represented by a real number, and

the one or more processors are configured to: in an updating process for each of the unit times, for each of the plurality of elements,

- update the first variable based on the second variable;
- when the first variable is smaller than a predetermined first value, change the first variable to the first value and change the second variable to a predetermined third value;
- when the first variable is greater than a predetermined second value greater than the first value, change the first variable to the second value and change the second variable to the third value; and
- add an acceleration value calculated by a predetermined computation to the second variable.

Technical Idea 2

The calculation device according to Technical Idea 1, wherein

the one or more processors are configured to: in the updating process for each of the unit times, for each of the plurality of elements,

- update the second variable before updating the first variable; and
- add the acceleration value to the second variable after updating the first variable.

Technical Idea 3

The calculation device according to Technical Idea 1, wherein

the one or more processors are configured to: in the updating process for each of the unit times, for each of the plurality of elements,

- update the second variable after updating the first variable; and
- add the acceleration value to the second variable after updating the first variable and before updating the second variable.

Technical Idea 4

The calculation device according to any one of Technical Ideas 1 to 3, wherein

the one or more processors are configured to:

- calculate, for each of the plurality of elements at the end time, a value of a discrete variable by binarizing the first variable with a preset threshold value, and
- output the calculated values of the discrete variables as a solution to the combinatorial optimization problem.

Technical Idea 5

The calculation device according to any one of Technical Ideas 1 to 4, wherein the acceleration value includes a value obtained by multiplying the corresponding second variable at a previous time by a product of a prescribed coefficient and the unit time.

Technical Idea 6

The calculation device according to Technical Idea 5, wherein the prescribed coefficient includes a predetermined constant.

Technical Idea 7

The calculation device according to Technical Idea 5, wherein the prescribed coefficient includes a value determined by a random number.

Technical Idea 8

The calculation device according to any one of Technical Ideas 1 to 7, wherein

the first value is −1,

the second value is +1, and

the third value is 0.

Technical Idea 9

The calculation device according to any one of Technical Ideas 1 to 8, wherein

each of the plurality of elements corresponds to one of a plurality of processing circuits, and

each of the processing circuits performs, for a corresponding element of the plurality of elements, a process of updating the first variable based on the second variable, processes of changing the first variable to the first value or the second value and changing the second variable to the third value, and a process of adding the acceleration value to the second variable.

Technical Idea 10

The calculation device according to any one of Technical Ideas 1 to 9, wherein

the one or more processors are configured to calculate, in the updating process for each of the unit times, for each of the plurality of elements, the first variable at a target time by adding a value obtained by multiplying the second variable, a predetermined constant, and the unit time to the first variable at a previous time by a unit time before the target time.

Technical Idea 11

The calculation device according to Technical Idea 10, wherein

the one or more processors are configured to calculate, in the updating process for each of the unit times, for each of the plurality of elements,

- calculate an external force, based on the first variable of each of the plurality of elements and a predetermined action coefficient for each of sets of a target element and a corresponding element of the plurality of elements; and
- calculate the second variable at the target time by adding a value obtained by multiplying a value determined by the external force and the first variable by the unit time, to the second variable at the previous time.

Technical Idea 12

The calculation device according to Technical Idea 11, wherein

the combinatorial optimization problem includes N discrete variables, and

the one or more processors are configured to calculate the first variable at the target time for an i^thelement corresponding to an i^thdiscrete variable of the N discrete variables, by Equation (101) or (102),

x_i(t+Δt)=x_i(t)+Dy_i(t)Δt (101)
x_i(t+Δt)=x_i(t)+Dy_i(t+Δt)Δt (102)

where

N is an integer equal to or greater than 2,

i is any integer of 1 to N,

D is the constant,

At is the unit time,

t is the previous time,

t+Δt is the target time,

x_i(t) is the first variable of the i^thelement at the previous time,

y_i(t) is the second variable of the i^thelement at the previous time,

x_i(t+Δt) is the first variable of the i^thelement at the target time, and

y_i(t+Δt) is the second variable of the i^thelement at the target time.

Technical Idea 13

The calculation device according to Technical Idea 12, wherein

the combinatorial optimization problem is a quadratic unconstrained binary optimization (QUBO) problem, and

the one or more processors are configured to calculate the second variable of the i^thelement at the target time by Equation (103) or Equation (104),

y_i(t+Δt)=y_i(t)+[{−D+p(t+Δt)}x_i(t+Δt)+f_i(t+Δt)]Δt (103)
y_i(t+Δt)=y_i(t)+[{−D+p(t)}x_i(t)+f_j(t)]Δt (104)

where f_i(t+Δt) is given by Equation (105), and f_i(t) is given by Equation (106),

f_i(t+Δt)=−cz_i(t+Δt) (105)
f_i(t)=−cz_i(t) (106)

where

z_i(t+Δt) is given by Equation (107), and

z_i(t) is given by Equation (108),

$\begin{matrix} z_{i} (t + Δ t) = - h_{i} α (t + Δ t) - \sum_{j = 1}^{N} J_{i, j} x_{j} (t + Δ t) & (107) \\ z_{i} (t) = - h_{i} α (t) - \sum_{j = 1}^{N} J_{i, j} x_{j} (t) & (108) \end{matrix}$

where

j is any integer of 1 to N,

h_iis an i^thlocal magnetic field coefficient included in a predetermined array including N local magnetic field coefficients,

J_i,jis a coupling coefficient at an i^throw and a j^thcolumn included in a predetermined matrix including N×N coupling coefficients,

c is a coefficient,

x_j(t) is the first variable at the previous time for a j^thelement corresponding to a j^thdiscrete variable of the N discrete variables,

x_j(t+Δt) is the first variable of the j^thelement at the target time,

p(t) is a predetermined function with t as a variable, in which p(t) increases as t increases, becomes 0 with t at the initial time, and becomes 1 with t at the end time, and

α(t) is a predetermined function with t as a variable.

Technical Idea 14

A calculation method of solving a combinatorial optimization problem by an information processing device, the calculation method including:

alternately updating, by the information processing device, for a plurality of elements each associated with a first variable and a second variable, the first variable and the second variable, sequentially for unit times from an initial time to an end time; and

outputting, by the information processing device, a solution to the combinatorial optimization problem based on the first variables of the plurality of elements at the end time, wherein

the plurality of elements correspond to a plurality of discrete variables of the combinatorial optimization problem,

the first variables and the second variables are each represented by a real number, and

in an updating process for each of the unit times by the information processing device, the updating includes: for each of the plurality of elements,

- updating the first variable based on the second variable;
- when the first variable is smaller than a predetermined first value, changing the first variable to the first value and changing the second variable to a predetermined third value;
- when the first variable is greater than a predetermined second value greater than the first value, changing the first variable to the second value and changing the second variable to the third value; and
- adding an acceleration value calculated by a predetermined computation to the second variable.

Technical Idea 15

A computer program product including a computer readable medium including programmed instructions for solving a combinatorial optimization problem by an information processing device, the instructions causing the information processing device to execute:

alternately updating, for a plurality of elements each associated with a first variable and a second variable, the first variable and the second variable, sequentially for unit times from an initial time to an end time; and

outputting a solution to the combinatorial optimization problem based on the first variables of the plurality of elements at the end time, wherein

the plurality of elements correspond to a plurality of discrete variables of the combinatorial optimization problem,

the first variables and the second variables are each represented by a real number, and

in an updating process for each of the unit times, the updating includes: for each of the plurality of elements,

- updating the first variable based on the second variable;
- when the first variable is smaller than a predetermined first value, changing the first variable to the first value and changing the second variable to a predetermined third value;
- when the first variable is greater than a predetermined second value greater than the first value, changing the first variable to the second value and changing the second variable to the third value; and
- adding an acceleration value calculated by a predetermined computation to the second variable.

Number	Name	Date	Kind
10860679	Goto et al.	Dec 2020	B2
20170104493	Goto	Apr 2017	A1
20190266212	Goto	Aug 2019	A1
20200090066	Kanao et al.	Mar 2020	A1
20210056160	Goto et al.	Feb 2021	A1

Number	Date	Country
2017-73106	Apr 2017	JP
2019-145010	Aug 2019	JP
2020-46784	Mar 2020	JP
2021-43667	Mar 2021	JP
WO 2021049063	Mar 2021	WO

Calculation device, calculation method, and computer program product

Information

Patent Number

Date Filed

Date Issued

Inventors

Original Assignees

Examiners

Agents

CPC

Field of Search

CPC

International Classifications

Term Extension

Abstract

Description

Claims

Priority Claims (1)

US Referenced Citations (5)

Foreign Referenced Citations (5)

Non-Patent Literature Citations (10)

Related Publications (1)

Entry
Zou, Yu, and Mingjie Lin. “Massively Simulating Adiabatic Bifurcations with FPGA to Solve Combinatorial Optimization: Proceedings of the 2020 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays.” ACM Conferences, Feb. 1, 2020, https://dl.acm.org/doi/10.1145/3373087.3375298. (Year: 2020).
A. Lucas, “Ising formulations of many NP problems,” Frontiers in Physics, vol. 2, Art. 5, 15 pages (2014).
M.W. Johnson et al., “Quantum annealing with manufactured spins,” Nature, vol. 473, pp. 194-198 (2011).
T. Inagaki et al., “A coherent Ising machine for 2000-node optimization problems,” Science, vol. 354, No. 6312, pp. 603-606 (2016).
H. Goto, “Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network,” Scientific Reports, 6:21686, DOI:10.1038/srep21686, 8 pages (2016).
M. Yamaoka et al., “A 20k-Spin Ising Chip to Solve Combinatorial Optimization Problems With CMOS Annealing,” IEEE J. Solid-State Circuits, vol. 51, No. 1, pp. 303-309 (2016).
H. Goto et al., “Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems,” Science Advances, 2019;5:eaav2372,8 pages (2019).
E.S. Tiunov et al., “Annealing by simulating the coherent Ising machine” Optics EXPRESS, vol. 27, No. 7, pp. 10288-10295 (2019).
H. Goto et al., “High-performance combinatorial optimization based on classical mechanics,” Science Advances, 2021;7:eabe7953, 9 pages (Feb. 3, 2021).
S. Tsukamoto et al., “An Accelerator Architecture for Combinatorial Optimization Problems,” Fujitsu Sci. Tech. J., vol. 53, No. 5, pp. 8-13 (2017).