JP7628866B2 - Information processing system, information processing method, and information processing program - Google Patents

Description

The present invention relates to an information processing system, an information processing method, and an information processing program.

Patent Document 1 discloses a semiconductor device having multiple units each including a first memory cell that stores a value expressing one spin of an Ising model in three or more states, a second memory cell that stores an interaction coefficient indicating an interaction from other spins that interacts with one spin, and a logic circuit that determines the next state of one spin based on a function having a value expressing the state of the other spins and the interaction coefficient as a constant or variable.

Patent document 2 also discloses a method for searching for an optimal solution for an Ising model with arbitrary connections by simultaneously probabilistically updating all spins while satisfying the theoretical background required by the Markov chain Monte Carlo method.

JP 2016-51314 A International Publication No. 2019/216277

Okuyama, T., Sonobe, T., Kawarabayashi, K. I., &Yamaoka, M. (2019). Binary optimization by momentum annealing. Physical ReviewE, 100(1), 012111. Botev, Z. I. (2017). The normal law under linearrestrictions: simulation and estimation via minimax tilting. Journal of theRoyal Statistical Society: Series B (Statistical Methodology), 79(1), 125-148. Neal, R. M. (1998). Suppressing random walks in Markovchain Monte Carlo using ordered overrelaxation. In Learning in graphical models(pp. 205-228). Springer, Dordrecht.

Many physical and social phenomena can be expressed by interaction models. An interaction model is defined by the multiple nodes that make up the model, the interactions between the nodes, and, if necessary, coefficients that act on each node. In the fields of physics and social science, various models, including the Ising model, have been proposed, and all of them can be interpreted as a form of interaction model.

Finding node states that minimize or maximize indicators related to this interaction model is important in solving social problems. Examples include the problem of detecting cliques in social networks and portfolio optimization problems in the financial sector. In the field of operations research, these are broadly classified as unconstrained binary quadratic programming problems and mixed binary quadratic programming problems.

The present invention has been made in consideration of the above-mentioned background, and aims to provide a technology that speeds up the search for optimal solutions to optimization problems, including ground state searches for interaction models.

In order to solve the above-mentioned problems, one aspect of the present invention is an information processing system that searches for an optimal solution to an optimization problem using a matrix product, and includes a processing unit that executes processing in cooperation with a storage unit, and the processing unit executes a conversion process that converts the optimization problem into a quadratic programming problem expressed by an objective function including a product of a first state variable vector, a predetermined matrix, and a second state variable vector, a state update process that performs probabilistic updates of the first state variable vector and the second state variable vector by repeating the processes of multiplying the predetermined matrix and the first state variable vector to calculate a first vector, performing probabilistic updates of the second state variable vector based on a probability distribution including a parameter based on the first vector, multiplying the predetermined matrix and the second state variable vector to calculate a second vector, and performing probabilistic updates of the first state variable vector based on a probability distribution including a parameter based on the second vector, and a function value calculation process that calculates the product by the inner product of the first state variable vector and the second vector to calculate a function value of the objective function.

According to the present invention, it is possible to speed up the search for an optimal solution to an optimization problem, including a ground state search for an interaction model.
Problems, configurations and effects other than those described above will become apparent from the following description of the preferred embodiment of the invention.

FIG. 1 is a conceptual diagram illustrating a relationship between a variable array and an objective function value of an optimization problem. FIG. 2 is a diagram illustrating an example of the relationship between variables of an objective function of an optimization problem. FIG. 2 is a diagram illustrating an example of the relationship between variables of an objective function of an optimization problem. FIG. 13 is a diagram illustrating an example of an optimal solution search process for multiple replicas. FIG. 2 is a diagram illustrating hardware of a computer that realizes the information processing device according to the embodiment. 1 is a functional block diagram illustrating an arithmetic circuit constituting an arithmetic device according to an embodiment; 1 is a functional block diagram illustrating a configuration of an information processing apparatus according to an embodiment. 11 is a flowchart illustrating an optimal solution search process executed by the information processing device according to the embodiment. 13 is a flowchart illustrating details of an interaction calculation process. FIG. 11 is a diagram illustrating a simulation result of a transition in the number of replicas due to resampling. 11A and 11B are diagrams illustrating simulation results comparing the calculation speed according to the embodiment with a comparative example.

The following describes an embodiment of the present invention with reference to the drawings. The embodiment is an example for explaining the present invention, and some parts have been omitted or simplified as appropriate for clarity of explanation. The present invention can also be implemented in various other forms. Unless otherwise specified, each component may be singular or plural.

When there are multiple components with the same or similar functions, they may be described using the same reference numerals with different subscripts. Also, when there is no need to distinguish between these multiple components, the subscripts may be omitted.

In some embodiments, the processing performed by executing a program is described. Here, the computer executes the program using a processor (e.g., CPU, GPU), and performs the processing defined by the program using storage resources (e.g., memory) and interface devices (e.g., communication ports). Therefore, the subject of the processing performed by executing the program may be the processor. Similarly, the subject of the processing performed by executing the program may be a controller, device, system, computer, or node having a processor. The subject of the processing performed by executing the program may be a calculation unit, and may include a dedicated circuit that performs specific processing. Here, the dedicated circuit is, for example, an FPGA (Field Programmable Gate Array), an ASIC (Application Specific Integrated Circuit), or a CPLD (Complex Programmable Logic Device).

The program may be installed on the computer from a program source. The program source may be, for example, a program distribution server or a computer-readable storage medium. When the program source is a program distribution server, the program distribution server may include a processor and a storage resource that stores the program to be distributed, and the processor of the program distribution server may distribute the program to be distributed to other computers. In addition, in an embodiment, two or more programs may be realized as one program, or one program may be realized as two or more programs.

(Theoretical Background of the Embodiments)
The optimization problem for optimizing an objective function associated with the interaction model has N variables s ₁ , ..., s _N , and the domain D _i of each variable s _i is either a binary value {-1, +1} or a continuous value [-1, +1]. Which the domain D _i of each variable s _i is is determined for each problem. The objective function H of the optimization problem is expressed by Equation 1. That is, the objective function H is expressed as a quadratic expression of the variable s.

式１において、ｓ＝［ｓ_１，…，ｓ_Ｎ］のＮ次元ベクトル、ＪはＮ×Ｎ対称行列、ｈはＮ次元ベクトルである。前述の通り、変数毎に定義域が異なるので、最適化問題は式２の通りに表される。

In formula 1, s=[s ₁ , . . . , s _N ] is an N-dimensional vector, J is an N×N symmetric matrix, and h is an N-dimensional vector. As described above, since the domain of definition differs for each variable, the optimization problem is expressed as in formula 2.

ここで、添字の集合Λ_ｂ、Λ_ｃを式３の通り定義する。

Here, the subscript sets Λ _b and Λ _c are defined as shown in Equation 3.

集合Ｓ_{ｍｉｘｅｄ}＝｛ｓ|ｓ_ｉ∈Ｄ_ｉ｝を定義する。これらの表記を用いると、式２は式４のように表現できる。

Define a set S _mixed ={s|s _i ∈D _i }. Using these notations, Equation 2 can be expressed as Equation 4.

以降、すべてのｉ∈Λ_ｂに対して行列Ｊのｉ行ｉ列目の要素は０とする。なぜならば、この変換は式２の最適解を変えないためである。

Hereinafter, for all _i∈Λb , the element in the i-th row and i-th column of the matrix J is set to 0, because this transformation does not change the optimal solution of Equation 2.

If D _i ={−1, +1} for all i, then this optimization problem is a combinatorial optimization problem called a ground state search problem of the Ising model. In this embodiment, in an optimization problem including a ground state search of the Ising model, an optimal solution or an approximate solution is searched for using an algorithm that utilizes Markov Chain Monte Carlo methods (hereinafter referred to as MCMC).

Figure 1 is a conceptual diagram illustrating the relationship between variable arrays and objective function values in an optimization problem, and shows the landscape of objective function values versus variable arrays. The horizontal axis of the graph is variable array s, and the vertical axis is the value of objective function H(s). MCMC repeats probabilistic transitions from the current state s to a state s' in the vicinity of state s. The probability of transitioning from state s to state s' is called the transition probability P(s, s'). Examples of transition probability P include the Metropolis method and the heat-bath algorithm.

The transition probability has a parameter called temperature, which indicates the ease of transition between states. When MCMC is performed while gradually decreasing the temperature from a large value, the objective function value asymptotically converges to the state with the lowest value (the state with the lowest energy, X in Figure 1). Methods that utilize this to find an optimal solution or an approximate solution to a minimization problem include simulated annealing (hereinafter referred to as SA (Simulated Annealing)) and momentum annealing (hereinafter referred to as MA (Momentum Annealing)), which was proposed in Non-Patent Document 1.

In solving the minimization problem shown in Equation 4, consider instead solving the minimization problem in Equation 5, where the set S _relaxed ={s|s _i ∈[-1, +1]}.

The optimal solution to equation 5 is expressed as s ^* = [ _s1 ^* , ..., _sN ^* ]. Although proof is omitted, s ⁺ = [ _s1 ⁺ , ..., _sN ⁺ ] obtained by equation 6 is one of the optimal solutions to equation 4. The goal of this embodiment is to search for the optimal solution to equation 2, but the desired solution s ⁺ can be obtained even if the transformation of equation 6 is obtained after the optimal solution s ^* of equation 5 is obtained. However, the function sgn is a function that returns +1 if the argument is 0 or more, and -1 otherwise.

Here, an N-dimensional vector v=[v _{1 ,} . . . , v _N ] is introduced to define a function H′ shown in Equation 7.

where function V(v) is defined as in Equation 8.

The matrix W=diag(w ₁ , . . . , w _N ) is an arbitrary diagonal matrix, and v _i is a real number ranging from [−1, +1]. Instead of the minimization problem shown in Equation 5, we introduce Equation 9, which is a minimization problem of H'(s, v).

Define two N-dimensional vectors x = s + v and y = s - v. The objective function of the optimization problem we want to solve is originally just H, but by introducing a function called V, we can obtain a new function that can be updated in parallel using MCMC. Then, the function H' can be rewritten as equation 10.

In other words, the minimization problem of Equation 5 can be rephrased as the minimization problem of Equation 11.

If the optimal solution of Equation 11 is expressed as x ^* and y ^* , then the following equality holds: s ^* = (x ^* + y ^* )/2. These arguments hold even if W is a zero matrix.

From the above, the optimal solution to the optimization problem expressed by Equation 2 can be found from the solution to the constrained quadratic programming problem shown in Equation 11. MCMC is used to find this solution.

FIG. 2 is a diagram illustrating the relationship between each variable of the objective function of the optimization problem, and shows a graphical model showing the relationship between each variable of the function G in Equation 11. FIG. 2 illustrates the case where N=6, but the same applies to general N. The relationship between each variable of the function G can be expressed by a complete bipartite graph. The variables that are multiplied by the variable x _i in the function G are only y ₁ , ..., y _N and x _i . When MCMC stochastically updates a variable value, it uses the value of the variable related to that variable. In other words, when updating the value of the variable x ₁ , y ₁ , ..., y _N and x ₁ are obtained, and other variables (here, x ₂ , ..., x _N ) are not referenced. This is also true for updating the value of other variables, for example, x ₂ . Therefore, if the value of the variable array y is constant, the theoretical requirements of MCMC are not violated even if the values of the array x are independently and simultaneously stochastically updated.

Similarly, the variables multiplied by the variable _yi are only _x1 , ..., _xN and _yi . Therefore, when the value of the variable array x is constant, each value of the array y can be independently and simultaneously stochastically updated.

As described above, by executing MCMC consisting of a procedure of repeating “simultaneous updating of _x1 , ..., _xN ” and “simultaneous updating of _y1 , ..., _yN ”, it is possible to search for arrays x and y that minimize the function G while enjoying the advantage of high speed achieved by parallelization.

Please note that in the discussion of this embodiment, no constraints are placed on matrix J. For example, even if all elements of matrix J are nonzero, the above discussion still holds, and parallel updates are possible.

Figure 3 is an example of a fully connected graph, illustrating the relationships between the variables of the objective function of an optimization problem. Figure 3 illustrates the case where N = 6, but the same is true for general values of N. On the other hand, when directly applying MCMC to the original problem, the minimization problem of Equation 2, the relationship of the variable array s is represented by a fully connected graph as shown in Figure 3, so the probability of only one variable at a time can be updated, and updates are limited to sequential updates.

From here, the procedure for probabilistic update for each variable will be described. The variable to be updated is denoted as x _i . When the values of variables y ₁ , . . . , y _N are constant, the existence probability p(x _i ) of variable x _i in the Boltzmann distribution of temperature T satisfies Equation 12.

Here, the variable A _i is the value calculated using equation 13.

Since variables x _i and y _i are |x _i |+|y _i |≦2, the range in which x _i can be moved is -(2-|y _i | ₎ or more and (2-|y _i | ₎ or less. Therefore, the next state of variable x _i can be sampled based on a truncated normal distribution with a domain of -(2-|y _i |) or more and (2-|y _i |) or less, in which the _normal distribution has a mean A i /w i and a variance T/w _i . This method means that the next state is determined regardless of the current state of x _i . The same applies to y _i . In this specification, when the variables x and y are not distinguished, they may be written as s.

Random numbers that follow the standard normal distribution can be generated using the Box-Muller method. Since the domain of definition is limited here, the algorithm shown in Non-Patent Document 2 can be used.

The search for an optimal solution can be regarded as sampling from an equilibrium state at temperature 0. Therefore, in order to realize a good solution search, it is preferable to converge to an equilibrium state in a short time. In order to improve the convergence to an equilibrium state, various techniques have been proposed in MCMC, and these can also be utilized. For example, Non-Patent Document 3 proposes an overrelaxation method. In this method, not only one state but K states are sampled from the Boltzmann distribution of temperature T as candidates for the next state. Then, in addition to the sampled K states, the total (K+1) states of the current state are rearranged and expressed as x _c ⁰ ≦...≦x _c ^r =x _i ≦x _c ^K. In other words, the current state is the (r+1)th value from the smallest of the (K+1) values. Then, x _c ^K+1-r is adopted as the next state. In this method, the next state depends on the current state of x _i .

As described above, the search for an optimal solution is performed by a stochastic search process that repeats stochastic transitions of the variables x and y that represent the state of the interaction model. However, as a result of the stochastic transitions, there are cases where a local solution is reached and an optimal solution or an approximate solution cannot be found even after many trials. Therefore, in this embodiment, multiple replicas corresponding to different states of the same model are created, and the stochastic search process is executed in parallel and independently for each replica. Then, the state is evaluated based on the objective function value, and resampling is performed to preferentially select and replicate replicas with high evaluations, thereby performing global state updates without falling into a local solution and speeding up the search for an optimal solution.

Figure 4 illustrates an example of an optimal solution search process for multiple replicas. In Figure 4, it is assumed that one replica is placed on each node.

For example, replicas 1-1, 1-2, and 1-3 arranged on three nodes #1, #2, and #3, each of which has its initial state set at temperature T=T ₀ (time step 0), are used to independently and in parallel search for an optimal solution. When temperature T=T _t (time step t), the objective function values (function G in Equation 11) of replicas 1-1, 1-2, and 1-3, which have been independently and in parallel search for an optimal solution up until that point, are calculated. The objective function of the j-th replica 1-j (j=1, ..., R) of the R replicas is defined as function G _j .

Next, the weight α _j of replica 1-j (j=1, . . . , R) is calculated as shown in Equation 14.

Here, Δβ _t is the difference (β _t -β _t-1 ) of the inverse temperature β _t (=1/T _t ) of the temperature T _t at the time step t when resampling is performed. From the weight α _j of this replica 1-j, the occurrence probability q _j representing the occurrence frequency of each replica 1-j after resampling is calculated as shown in Equation 15.

However, the denominator on the right side of Equation 15 represents the sum of the weights α _j of all replicas 1-j before resampling, that is, Equation 15 normalizes the weights α _j .

Then, the selection number _Nj of each replica 1-j is determined as the replica after resampling as shown in Equation 16. The selection number _Nj is the number of replicas 1-j used in the optimum solution search of the next time step. The decimal part of _Nj is rounded off appropriately. R represents the total number of replicas, which is constant before and after resampling, but may increase or decrease after resampling.

Then, the copy destination of the replica after resampling is determined according to the selection number _Nj , and the information of each replica 1-j is copied as necessary.

In the example of FIG. 4, R=3. Replica 1-1 is R×q ₁ =0.6, replica 1-2 is R×q ₂ =2.1, and replica 1-3 is R×q ₃ =0.03, so by rounding, one replica 1-1, two replicas 1-2, and zero replicas 1-3 are selected. That is, the replicas after resampling are one replica 1-1 and two replicas 1-2. Replica 1-1 and the first replica 1-2 are used continuously from before resampling. And, since the second replica 1-2 cannot be placed on node #2, the information of replica 1-2 is transferred to node #3 where replica 1-3 is not selected. After resampling, a stochastic search process is executed in parallel and independently for each of one replica 1-1 and two replicas 1-2.

When searching for an optimal solution, the replicas are executed in parallel, and when an execution termination condition is met, such as the temperature dropping to a certain level, an optimal solution or an approximate solution is output based on a variable representing the state of the replica with the smallest value of function G among the multiple replicas.

Based on the above, the configuration of the information processing device 10 according to this embodiment will be described with reference to Figures 5 to 7.

(Hardware configuration of information processing device 10)
5 is a diagram illustrating hardware of a computer that realizes an information processing device 10 according to an embodiment. The information processing device 10 has, as hardware, a processor 11, a main storage device 12, an auxiliary storage device 13, an input device 14, an output device 15, a communication device 16, one or more arithmetic devices 17, and a system bus 18 that communicatively connects these devices. The information processing device 10 may be realized, for example, using a virtual information processing resource such as a cloud server provided in part or in whole by a cloud system. The information processing device 10 may also be realized, for example, by a plurality of information processing devices that are communicatively connected and operate in cooperation with each other.

The processor 11 is configured using, for example, a CPU (Central Processing Unit) and an MPU (Micro Processing Unit). The main memory device 12 is a device that stores programs and data. For example, it is a ROM (Read Only Memory), SRAM (Static Random Access Memory), NVRAM (Non Volatile RAM), mask ROM (Mask Read Only Memory), PROM (Programmable ROM), RAM (Random Access Memory), DRAM (Dynamic Random Access Memory), etc. The auxiliary memory device 13 is a HDD (Hard Disk Drive), flash memory, SSD (Solid State Drive), optical memory device (CD (Compact Disc), DVD (Digital Versatile Disc), etc.), etc. The programs and data stored in the auxiliary memory device 13 are loaded into the main memory device 12 as needed.

The input device 14 is a user interface that accepts information input from the user, and is, for example, a keyboard, a mouse, a card reader, a touch panel, etc. The output device 15 is a user interface that provides information to the user, and is, for example, a display device (LCD (Liquid Crystal Display), a graphics card, etc.) that visualizes various information, an audio output device (speaker), a printer, etc. The communication device 16 is a communication interface that communicates with other devices, and is, for example, a NIC (Network Interface Card), a wireless communication module, a USB (Universal Serial Interface) module, a serial communication module, etc.

The arithmetic unit 17 is a device that executes processing related to searching for optimal solutions to combinatorial optimization problems. The arithmetic unit 17 may take the form of an expansion card that is attached to the information processing device 10, such as a GPU (Graphics Processing Unit). The arithmetic unit 17 is configured with hardware, such as a CMOS (Complementary Metal Oxide Semiconductor) circuit, FPGA, or ASIC.

The arithmetic unit 17 includes a control device, a storage device, an interface for connecting to the system bus 18, and transmits and receives commands and information to and from the processor 11 via the system bus 18. The arithmetic unit 17 may be communicatively connected to other arithmetic units 17 via a communication line, for example, and may operate in cooperation with the other arithmetic units 17. The functions realized by the arithmetic unit 17 may be realized, for example, by having a processor (CPU, GPU, etc.) execute a program.

In the information processing device 10, a program is read from the auxiliary storage device 13 and executed by the processor 11 and the main storage device 12 in cooperation with each other, thereby realizing the functions of the information processing device 10. Alternatively, the program for realizing the functions of the information processing device 10 may be obtained from an external computer by communication via the communication device 16. Alternatively, the program for realizing the functions of the information processing device 10 may be recorded on a portable recording medium (optical disk, magnetic disk, magneto-optical disk, semiconductor storage medium, etc.) and read by a medium reading device.

In addition, the information processing device 10 may be an information processing system in which multiple devices work together to execute processing. A program that realizes the information processing system realizes functions similar to those of the information processing device 10 by having each device realize each function through each program.

(Configuration of arithmetic circuit 500)
FIG. 6 is a functional block diagram illustrating an arithmetic circuit 500 constituting the arithmetic device 17 according to the embodiment, and shows the operation principle of the arithmetic device 17. Each arithmetic device 17 is configured to include one or more arithmetic circuits 500. Each arithmetic circuit 500 corresponds to one replica. In this embodiment, the arithmetic circuit 500 is implemented by software using a number of arithmetic units mounted on a GPU, vector type computer, etc., for example, a difference calculation block 514, a sampling block 515, a next state determination block 516, etc. However, this is not limited to this, and may be implemented by hardware. The arithmetic circuit 500 realizes a function of sampling a variable array x ₁ , ..., x _N or a variable array y ₁ , ..., y _N from a Boltzmann distribution (Equation 12) at temperature T.

The arithmetic circuit 500 includes a variable memory 511, a nonlinear coefficient memory 512, a linear coefficient memory 513, a difference calculation block 514, a sampling block 515, and a next state determination block 516.

The variable memory 511 of each arithmetic circuit 500 stores information indicating the above-mentioned variables x ₁ , . . . , _xN and y1, _.

The nonlinear coefficient memory 512 stores information representing the matrix J. The matrix J is generally a symmetric matrix, and this symmetry can be used to reduce the amount of usage of the nonlinear coefficient memory 512. The linear coefficient memory 513 stores information representing the vector h.

The control signal EN, the weight signal SW, and the temperature signal TE are input to the calculation device 17.

The control signal EN is a signal that periodically repeats the values of H (high) and L (low), and indicates whether the variable array x or y is being updated. For example, when EN is H, the variable array x is updated, and when EN is L, y is updated. This control signal EN simultaneously updates variables _x1 , ..., _xN , and simultaneously updates variables _y1 , ..., _yN . For simplification in Fig. 5, the control signal EN is input only to sampling block 515, but it is also input in the same way to other locations that require this signal, such as variable memory 511.

The weight signal SW is a signal that represents a vector of N elements that represent the diagonal components of the diagonal matrix W.

The difference calculation block 514 receives as input the value of matrix J stored in nonlinear coefficient memory 512, vector h stored in linear coefficient memory 513, weight signal SW, and variable s (x or y) stored in variable memory 511. The difference calculation block 514 outputs y+h (J+diag( _w1 , ..., _wN )) when the control signal EN is H, and outputs x+h (J+diag( _w1 , ..., _wN )) when EN is L. This output value corresponds to the above-mentioned _Ai .

Sampling block 515 receives the output of difference calculation block 514, weight signal SW, temperature signal TE that holds the value of the temperature parameter, control signal EN, and values of other variables. Then, as the i-th element, it randomly samples and outputs from a truncated normal distribution expressed by Equation 12, whose domain is from -(2-|y _i |) to (2-|y _i |) when control signal EN is H, and from -(2-|x _i |) to (2-|x _i |) when EN is L.

The next state determination block 516 determines the next state of the variable based on one or more values output from the sampling block 515. If the MCMC update rule is set to a simple heat bath method, the next state determination block 516 only needs to receive one output value from the sampling block 515 and write it directly to the variable memory 511. If the well-known over-relaxation method is used as the MCMC update rule, the next state determination block 516 receives multiple values from the sampling block 515 and the current value of the variable to be updated from the variable memory 511, selects one according to the over-relaxation method, and writes it to the variable memory 511. As is well known, in the over-relaxation method, the next state is determined so that the correlation with the previous state is negative.

(Functional configuration of information processing device 10)
7 is a functional block diagram illustrating the configuration of the information processing device 10 according to the embodiment, and shows main functions (software configuration) of the information processing device 10. The information processing device 10 includes a processing unit 101, a storage unit 102, and one or more arithmetic units 17.

The processing unit 101 has a model conversion unit 101a, a model coefficient setting unit 101b, a weight setting unit 101c, a variable value initialization unit 101d, a temperature setting unit 101e, an interaction calculation execution unit 101f, and a variable value reading unit 101g. These functions are realized by the processor 11 working in cooperation with the main memory unit 12 to read and execute a program stored in the auxiliary memory unit 13, or by hardware provided in the calculation unit 17. In addition to the above functions, the information processing device 10 may also have other functions such as an operating system, a file system, a device driver, a DBMS (Data-Base Management System), etc.

The memory unit 102 is realized by the main memory device 12 or the auxiliary memory device 13, and stores problem data 102a, quadratic programming format problem data 102b, domain data 102c, and a calculation device control program 102d. The problem data 102a is, for example, data in which an optimization problem or the like is described in a known, predetermined description format. The problem data 102a is, for example, set by a user via a user interface (input device, output device, communication device, etc.).

The quadratic programming problem data 102b is data generated by the model conversion unit 101a by converting the problem data 102a into data in a format that matches the format of the quadratic programming problem shown in Equation 4. In this conversion, the domain of each given variable is written to the domain data 102c. The domain data 102c indicates, for example, whether each variable takes on a binary value or a real number value. The calculation device control program 102d is a program that is executed when the interaction calculation execution unit 101f controls the calculation device 17, or that is loaded by the interaction calculation execution unit 101f into each calculation device 17 and executed by the calculation device 17.

The model conversion unit 101a converts the problem data 102a into quadratic programming problem data 102b, which is a format for quadratic programming problems. For this purpose, a function for deriving equation 11 from equation 1 may be implemented in the model conversion unit 101a as software or hardware. The function of the model conversion unit 101a does not necessarily have to be implemented in the information processing device 10, and the information processing device 10 may import quadratic programming problem data 102b generated by another information processing device or the like via the input device 14 or communication device 16.

The model coefficient setting unit 101b, weight setting unit 101c, variable value initialization unit 101d, and temperature setting unit 101e generate R calculation circuits 500 in the calculation device 17. Then, the following process is performed on the R calculation circuits 500 to create R replicas 1-j (j=1, 2, ..., R). The R calculation circuits 500 are placed in one or more calculation devices 17.

The model coefficient setting unit 101b sets the matrix J of Equation 4 in the nonlinear coefficient memory 512 and the vector h in the linear coefficient memory 513 based on the quadratic programming problem data 102b.

The variable value initialization unit 101d initializes the values of each variable stored in the variable memory 511 of the arithmetic circuit 500. The variable value initialization unit 101d may, for example, uniformly determine the value of each variable by random sampling from -1 to +1. At this time, care must be taken to satisfy the constraint on the variables, |x _i |+|y _i |≦2. It should also be noted that the values of each variable at this time are treated as continuous values.

The temperature setting unit 101e sets the temperature T that the interaction calculation execution unit 101f uses when searching for an optimal solution.

The interaction calculation execution unit 101f causes each calculation circuit 500 of each replica 1-j to execute a calculation (hereinafter referred to as an interaction calculation process) to search for variable arrays x and y that minimize the function G expressed by equation 11 for each temperature T set by the temperature setting unit 101e. During the interaction calculation process, the interaction calculation execution unit 101f changes the temperature T, for example, from a higher temperature to a lower temperature.

Specifically, the probabilistic update of the variable array in each replica 1-j uses variable y when updating variable x, and variable x when updating variable y, as shown in equations 12 and 13. Specifically, as can be seen from equations 12 and 13, the multiplication result of matrix J and variable y is used when updating variable x. Similarly, the multiplication result of matrix J and variable x is used when updating variable y.

The interaction calculation execution unit 101f also executes resampling at a predetermined period. That is, the interaction calculation execution unit 101f calculates the energy (objective function value (value of function G _j )) for each replica 1-j created in each calculation circuit 500 (node) at temperature T of time step t.

Here, the multiplication of the matrix J of the function G _j by the variable y (the first term on the right side of Equation 11) has a large calculation load. However, this multiplication of the matrix J by the variable y has been performed in the previous update step of the variable x, and the multiplication result Jy is stored in the work area. Therefore, when calculating the value of the function G _j , instead of multiplying the matrix J by the variable y, the multiplication result Jy stored in the work memory is used. Since the multiplication result Jy is a vector, the calculation of the first term on the right side of the function G _j is replaced with an inner product calculation with the vector x ^T. In this way, by using the calculation results of the update process of the variables x and y to calculate the objective function value, the calculation amount of the objective function value can be suppressed, and the calculation speed of the optimization search can be increased.

Then, the interaction calculation execution unit 101f determines the number N _j of selected replicas 1-j after resampling using Equations 14, 15, and 16 based on the energy of each replica 1-j at time step t.

Then, the interaction calculation execution unit 101f determines the copy destination of each replica 1-j of the increment after resampling among the selected number N _j of each replica 1-j. The copy destination of each replica 1-j is determined preferentially to the same calculation device 17 (own node) as the copy source replica 1-j within the upper limit of the own node. If the upper limit of the own node is exceeded, the copy destination of the replica is randomly determined to another calculation device 17 (other node) different from the own node within each upper limit. The upper limit is the resource upper limit of each node, or the number of replicas arranged in each node before resampling. Then, the information of the copy source replica is transferred to the copy destination.

In this way, by giving priority to the copy destination within the own node and avoiding copying to other nodes where communication occurs when copying replica information, it is possible to reduce waiting times for information transfer and prevent a decrease in calculation speed.

When the interaction calculation execution unit 101f finishes searching for the optimal solution, the variable value reading unit 101g reads out the variable arrays x and y stored in the variable memory 511 of the calculation circuit 500 corresponding to the replica with the smallest value of the function G. The value read out here is the solution of Equation 11. According to the above discussion, the N-dimensional vector s*=(x+y)/2 is calculated. Then, the domain data 102c is read out, and the vector s+ obtained by Equation 6 is output to the output device 15 or the communication device 16 as the final solution. In other words, if the i-th domain in the domain data 102c is found to be {-1, +1}, sgn(s*i) is output, and if the i-th domain is [-1, +1], s _i itself is output. In this way, a solution according to the defined value range is found.

(Flow of optimal solution search process)
8 is a flowchart illustrating an optimum solution search process executed by the information processing apparatus 10 according to the embodiment. The optimum solution search process is started by receiving an instruction from a user via the input device 14, for example.

In S11, the model conversion unit 101a converts the problem data 102a into quadratic programming problem data 102b. In the quadratic programming problem data 102b, for example, the matrix J and vector h in the objective function H expressed in Equation 1 are expressed in an arbitrary format. If the storage unit 102 has already stored the quadratic programming problem data 102b, S11 is omitted. The process of S11 and the process from S12 onwards may be executed by different devices. Furthermore, the process of S11 and the process from S12 onwards may be executed at different times (for example, the process of S11 may be executed in advance).

The processes of S12 to S15 are performed for R replicas created by R arithmetic circuits 500 of the arithmetic device 17.

In S12, the model coefficient setting unit 101b sets the values of the matrix J and the vector h in the nonlinear coefficient memory 512 and the linear coefficient memory 513. The memory values can be set or edited by the user via a user interface (realized, for example, by the input device 14, the output device 15, the communication device 16, etc.).

In S13, the weight setting unit 101c determines the value of the weight signal SW. As explained above in Equation 8, the weight signal SW is allowed to take any value when searching for an optimal solution. Therefore, the signal value may always be 0. In this case, the calculation load can be reduced. Also, as shown in Equations 3 to 5 of Patent Document 2, the weight signal SW may be determined from the eigenvalues of matrix J. Alternatively, the weight signal SW may be determined from the row sum of matrix J. The calculation for calculating the value of the weight signal SW may be performed in the arithmetic unit 17 or the processor 11. Alternatively, the user may set it himself (S13).

In S14, the variable value initialization unit 101d initializes the values of each variable stored in the variable memory 511. The values stored in the variable memory 511 are continuous values. As mentioned above, the initial values may be random. With the above, the parameters expressing Equation 11 have been set.

In S15, the temperature setting unit 101e initializes a series T _γ (γ=1, 2, 3, ...) of temperature parameters to be used in the optimum solution search. Note that the above subscript γ indicates the type of the set temperature T. The method of setting the temperature T can be, for example, the method described in Patent Document 1.

In S16, the interaction calculation execution unit 101f executes the interaction calculation process. Details of the interaction calculation process will be described later with reference to FIG. 9.

In S17, the interaction calculation execution unit 101f determines whether the stop condition is met (for example, whether the temperature T has reached a preset minimum temperature). If the interaction calculation execution unit 101f determines that the stop condition is met (S17: Yes), it proceeds to S18, whereas if it determines that the stop condition is not met (S17: No), it returns to S16.

In S18, the variable value reading unit 101g reads the values of the variables stored in the variable memory 511 from the arithmetic circuit 500 of the replica 1-j having the smallest value of the function G _j . The variable value reading unit 101g also reads the domain of each variable of the quadratic programming problem data 102b stored in the domain data 102c. The variable value reading unit 101g then calculates a vector through a transformation based on Equation 6, and outputs it as a solution to Equation 2 or Equation 4. This completes the optimal solution search process.

(Detailed flow of interaction calculation process)
Fig. 9 is a flowchart illustrating details of the interaction calculation process in S16 in Fig. 8. In S161, the interaction calculation execution unit 101f executes probabilistic simultaneous updating of the variable arrays of each replica 1-j through calculations by each calculation circuit 500.

In S162, the interaction calculation execution unit 101f calculates the value of the function G _j (Equation 11) for each replica 1-j. As described above, when calculating the objective function value of the function G _j for each replica 1-j, the calculation results of the update process of the variables x and y are used. Note that the calculation of the value of the function G _j may be performed by the calculation circuit 500 of each replica 1-j.

In S163, the interaction calculation execution unit 101f determines the selection number _Nj of each replica 1-j and the copy destination of the increased replica.

In S164, the interaction calculation execution unit 101f transmits information about the replica to be copied to the destination node (the node itself or another node) determined in S163.

(Simulation results of the change in the number of replicas)
FIG. 10 is a diagram illustrating a simulation result of the transition of the number of replicas due to resampling. In FIG. 10, the vertical axis represents time, and the horizontal axis represents the number of replicas. FIG. 10 shows how replicas 1-1, 1-2, 1-3, 1-4, 1-5, 1-6, 1-7, 1-8, 1-9, and 1-10, each of which was one at time 0 with the highest temperature, increase and decrease through repeated resampling. In FIG. 10, replicas with the same original replica at time 0 are all in the same system. Therefore, it should be noted that replicas with the same filled or hatched pattern after time 1 do not necessarily have the same variables. For example, replicas with some excellent objective function values, such as replicas 1-6 and 1-9 at time 50 with the lowest temperature, are predominant, which increases the possibility that the objective function value will converge to the lowest state more quickly, and the optimization problem can be processed more quickly.

(Effects of the Disclosed Technology Including the Embodiments)
In the disclosed technology including the above-mentioned embodiment, the optimization problem is converted into a quadratic programming problem represented by an objective function including a product of a first state variable vector, a predetermined matrix, and a second state variable vector, and when calculating the function value of the objective function, the objective function is calculated by the inner product of the first state variable vector and the second vector. Here, the second vector is a second vector obtained by performing the stochastic update of the first state variable vector and the second state variable vector by repeating the process of multiplying the first state variable vector by a predetermined matrix and the first state variable vector, performing the stochastic update of the second state variable vector based on a probability distribution including a parameter based on the first vector, multiplying the second state variable vector by a predetermined matrix and the second state variable vector, and performing the stochastic update of the first state variable vector based on a probability distribution including a parameter based on the second vector, and the second vector is stored in the storage unit. Therefore, by reusing the result of the calculation process that appears in the iterative calculation of the optimal solution search, the objective function value can be calculated at low cost, and the process of the optimal solution search can be accelerated.

In addition, in the disclosed technology, state updates and function value calculations are performed independently and in parallel for multiple replicas corresponding to the first state variable vector and the second state variable vector, so state updates can be performed globally without being limited to individual local solutions.

In addition, in the disclosed technology, state updates and function value calculations for each replica are performed by multiple arithmetic units, so that parallel processing can be performed while distributing the load, thereby speeding up processing.

In addition, in the disclosed technology, a function value for each replica is calculated by a function value calculation process, a weight for each replica is calculated based on the difference between the function value for each replica and the inverse temperature, and resampling is performed to determine the number of replicas to select from the multiple replicas based on the weight. Therefore, by reusing the results of the calculation process that appears in the iterative calculation of the optimal solution search, the objective function value can be calculated at a low cost on the order of O( _log2n ), which is a heavy calculation for calculating the function value each time resampling is repeated, and the processing of the optimal solution search can be accelerated.

In addition, with the disclosed technology, the number of replicas selected in resampling is the same as the total number of replicas before resampling, so limited resources can be allocated more to the series of replicas with higher objective function values, efficiently speeding up the search for the optimal solution.

In addition, in the disclosed technology, the number of replicas selected in resampling is greater than the total number of replicas before resampling, so the optimum solution search can be performed using many system replicas using abundant resources, further speeding up the optimum solution search.

In addition, with the disclosed technology, the number of replicas selected in resampling is smaller than the total number of replicas before resampling, so that limited resources can be allocated to the series of replicas with high objective function values, efficiently speeding up the search for an optimal solution.

Furthermore, in the disclosed technology, depending on the number of replicas selected in resampling, if a copy of the replica can be placed on the arithmetic device on which the selected replica is placed, this arithmetic device is determined as the copy destination for placing the copy of the replica preferentially. On the other hand, if placement is not possible, another arithmetic device is determined as the copy destination. This reduces the occurrence of inter-node communication when transferring (copying) replica information during resampling, thereby reducing communication time and speeding up the search for an optimal solution.

As examples of the above-mentioned embodiment, the simulation results of the calculation speed for Example 1: implementation with 1 GPU (implementation of one calculation unit 17 (GPU)) and Example 2: implementation with 4 GPUs (implementation of four calculation units 17 (GPU)) are shown below.

Figure 11 is a diagram illustrating simulation results comparing the calculation speed according to the embodiment with a comparative example. In Figure 11, the horizontal axis represents the calculation time for one annealing run, and the vertical axis represents the objective function value. A total of 100 trials were run for each run with different random number seeds, and the average, first quartile, and third quartile of the objective function value at each time were plotted. The dashed line indicates the average value, the lower outline of the gray line that overlaps with the dashed line indicates the first quartile, and the upper outline indicates the third quartile.

As shown in Table 1, in Comparative Example 1, SA (Simulated Annealing) with 50 nodes was performed twice for a total of 100 trials. In Comparative Example 2, MA (Momentum Annealing) without PMC (Parallel Markov Chains) was performed 100 times with 1 GPU. In Comparative Example 3, MA without PMC was performed 100 times with 4 GPU. In Example 1, MA with PMC was performed 100 times with 1 GPU. In Example 2, MA with PMC was performed 100 times with 4 GPU.

11, Examples 1 and 2 reached near the minimum energy state (objective function value = -4.40 x 10 ⁵ ) more quickly compared to any of Comparative Examples 1 to 3. Moreover, Example 1 reached near the minimum energy state more quickly compared to Comparative Example 2, which also had 1 GPU implementation. Moreover, Example 2 reached near the minimum energy state more quickly compared to Comparative Example 3, which also had 4 GPU implementation. Moreover, Example 2 reached near the minimum energy state more quickly compared to Example 1.

The above simulation results show that, even with the same MA, a 4-GPU implementation is more advantageous in terms of calculation speed than a 1-GPU implementation, and furthermore, in the MA with PMC of this embodiment, the more GPUs implemented, the more advantageous it is in terms of calculation speed.

Other Embodiments
The arithmetic circuit 500 may be configured with software or hardware as long as it has a function of executing calculations to solve the optimization problem described above. Specifically, not only the hardware implemented with electronic circuits (digital circuits, etc.) in the annealing method, but also the method implemented with superconducting circuits, etc. may be used. Also, hardware that realizes the Ising model by a method other than the annealing method may be used. For example, a laser network method (optical parametric oscillation), a quantum neural network, etc. are known. In addition, although the concept is partially different, a quantum gate method in which the calculations performed in the Ising model are replaced with gates such as a Hadamard gate, a rotation gate, and a controlled NOT gate can also be adopted as a configuration of this embodiment.

As an example of implementation of the arithmetic circuit 500, it can be implemented as a CMOS (Complementary Metal-Oxide Semiconductor) integrated circuit or a logic circuit on an FPGA. For example, as disclosed in Patent Document 1, a large number of units that apply SRAM (Static Random Access Memory) technology may be arranged, and each unit may be provided with a memory for storing variables and a circuit for updating the variables.

The present invention is not limited to the above-described embodiments, but includes various modifications. For example, the above-described embodiments have been described in detail to clearly explain the present invention, and are not necessarily limited to those having all of the configurations described. Furthermore, it is possible to replace part of the configuration of one embodiment with the configuration of another embodiment, and to add the configuration of another embodiment to the configuration of one embodiment, so long as there is no contradiction. It is also possible to add, delete, replace, integrate, or distribute part of the configuration of each embodiment. Furthermore, the configurations and processes shown in the embodiments can be appropriately distributed, integrated, or replaced based on processing efficiency or implementation efficiency.

10: Information processing device, 11: Processor, 12: Main storage device, 13: Auxiliary storage device, 17: Arithmetic device, 511: Variable memory, 512: Nonlinear coefficient memory, 513: Linear coefficient memory, 514: Difference calculation block, 515: Sampling block, 516: Next state determination block, 101: Processing unit, 101a: Model conversion unit, 101b: Model coefficient setting unit, 101c: Weight setting unit, 101d: Variable value initialization unit, 101e: Temperature setting unit, 101f: Interaction calculation execution unit, 101g: Variable value reading unit, 102: Memory unit, 102a: Problem data, 102b: Quadratic programming format problem data, 102c: Domain data, 102d: Arithmetic device control program

Claims (14)

An information processing system for searching for an optimal solution to an optimization problem using matrix multiplication, comprising:
A processing unit that cooperates with the storage unit to execute processing,
The processing unit includes:
a conversion process for converting the optimization problem into a quadratic programming problem expressed by an objective function including a product of a first state variable vector, a predetermined matrix, and a second state variable vector;
a state update process for performing probabilistic updates of the first state variable vector and the second state variable vector by repeating a process of multiplying the first state variable vector by the predetermined matrix to calculate a first vector, performing probabilistic updates of the second state variable vector based on a probability distribution including a parameter based on the first vector, multiplying the second state variable vector by the predetermined matrix to calculate a second vector, and performing probabilistic updates of the first state variable vector based on a probability distribution including a parameter based on the second vector;
a function value calculation process for calculating a function value of the objective function by calculating the product by an inner product of the first state variable vector and the second vector. 2. The information processing system according to claim 1,
The processing unit includes:
an information processing system comprising: an information processing unit that executes the state update process and the function value calculation process independently and in parallel for each of a plurality of replicas corresponding to the first state variable vector and the second state variable vector; 3. The information processing system according to claim 2,
A plurality of computing devices;
The processing unit includes:
the state update process and the function value calculation process for each replica are executed by the plurality of arithmetic units. 4. The information processing system according to claim 3,
The processing unit includes:
calculating the function value for each replica by the function value calculation process;
a weight for each replica is calculated based on the difference between the function value and the inverse temperature for each replica, and a resampling process is performed to determine the number of replicas to be selected from the plurality of replicas based on the weight. 5. The information processing system according to claim 4,
The information processing system according to claim 1, wherein the selected number is equal to a total number of the plurality of replicas. 5. The information processing system according to claim 4,
The information processing system, wherein the selected number is greater than a total number of the plurality of replicas. 5. The information processing system according to claim 4,
The information processing system according to claim 1, wherein the selected number is smaller than a total number of the plurality of replicas. 5. The information processing system according to claim 4,
The processing unit includes:
an information processing system, characterized in that in the resampling process, depending on the selection number, if a copy of the replica can be placed on the arithmetic device on which the selected replica is placed, the arithmetic device is determined as the destination for placing the copy of the replica, and if placement is not possible, another of the arithmetic devices is determined as the destination. An information processing method performed by an information processing system that searches for an optimal solution to an optimization problem using matrix multiplication, comprising the steps of:
a conversion process for converting the optimization problem into a quadratic programming problem expressed by an objective function including a product of a first state variable vector, a predetermined matrix, and a second state variable vector;
a state update process for performing probabilistic updates of the first state variable vector and the second state variable vector by repeating a process of multiplying the first state variable vector by the predetermined matrix to calculate a first vector, performing probabilistic updates of the second state variable vector based on a probability distribution including a parameter based on the first vector, multiplying the second state variable vector by the predetermined matrix to calculate a second vector, and performing probabilistic updates of the first state variable vector based on a probability distribution including a parameter based on the second vector;
a function value calculation process of calculating a function value of the objective function by calculating the product by an inner product of the first state variable vector and the second vector. 10. The information processing method according to claim 9,
The state update process and the function value calculation process include:
an information processing method, characterized in that the information processing method is executed independently and in parallel for each of a plurality of replicas corresponding to the first state variable vector and the second state variable vector. 11. The information processing method according to claim 10,
The information processing system includes a plurality of computing devices,
the state update process and the function value calculation process for each replica are executed by the plurality of arithmetic units. The information processing method according to claim 11,
calculating the function value for each replica by the function value calculation process;
a resampling process for calculating a weight for each replica based on the difference between the function value and the inverse temperature for each replica, and determining the number of replicas to be selected from the plurality of replicas based on the weight. 13. The information processing method according to claim 12,
an information processing method, characterized in that in the resampling process, depending on the selection number, if a copy of the replica can be placed on the arithmetic device on which the selected replica is placed, the arithmetic device is determined as the copy destination for placing the copy of the replica, and if placement is not possible, another of the arithmetic devices is determined as the copy destination. A program for causing a computer to function as the information processing system according to any one of claims 1 to 8.

Images