Information processing system, information processing method, and information processing program

Abstract

To provide an information processing system that speeds up an optimum solution search for an optimization problem including a ground state search for an interaction model.SOLUTION: An information processing system for searching for an optimal solution for an optimization problem using a matrix product, is configured to: stochastically update a second state variable vector on the basis of a probability distribution which contains a parameter based on a first vector obtained by multiplying a predetermined matrix by a first variable vector, of an object function which contains a product of the first state variable vector of a quadratic programming problem converted from an optimization problem, the predetermined matrix, and the second state variable vector; repeat processing of stochastically updating the first state variable vector on the basis of a probability distribution which contains a parameter based on a second vector obtained by multiplying the predetermined matrix by the second state variable vector; and stochastically update the first and second state variable vectors. The system calculates a function value of the object function by calculating a product by an inner product of the first state variable vector and the second vector.SELECTED DRAWING: Figure 9

Description

[Technical field]

[0001] The present invention relates to an information processing system, an information processing method, and an information processing program. [Background technology]

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

[0003] Furthermore, Patent Document 2 discloses a method for searching for an optimal solution for an Ising model having any coupling by simultaneously probabilistically updating all spins while satisfying the theoretical background required by the Markov chain Monte Carlo method. [Prior art documents] [Patent documents]

[0004] [Patent Document 1] JP 2016-51314 A [Patent Document 2] International Publication No. 2019 / 216277 [Non-patent literature]

[0005] [Non-Patent Document 1] Okuyama, T., Sonobe, T., Kawarabayashi, KI, &Yamaoka, M. (2019). Binary optimization by momentum annealing. Physical ReviewE, 100(1), 012111. [Non-Patent Document 2] Botev, ZI (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. [Non-Patent Document 3] Neal, RM (1998). Suppressing random walks in Markovchain Monte Carlo using ordered overrelaxation. In Learning in graphical models(pp. 205-228). Springer, Dordrecht. Summary of the Invention [Problem to be solved by the invention]

[0006] Many physical and social phenomena can be expressed by interaction models. An interaction model is defined by multiple nodes that make up the model, the interactions between the nodes, and, if necessary, coefficients acting 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.

[0007] Finding node states that minimize or maximize indicators related to this interaction model is important in solving social problems. For example, problems such as detecting cliques in social networks and portfolio optimization problems in the financial sector can be categorized as unconstrained binary quadratic programming problems and mixed binary quadratic programming problems in the field of operations research.

[0008] The present invention has been made in view of the above-mentioned background, and has an object to provide a technique for speeding up an optimal solution search for an optimization problem, including a ground state search for an interaction model. [Means for solving the problem]

[0009] 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, the information processing system having a processing unit that executes processing in cooperation with a storage unit, the processing unit executes the following processes: 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: calculating a first vector by multiplying the predetermined matrix by the first state variable vector, performing probabilistic updating of the second state variable vector based on a probability distribution including a parameter based on the first vector, multiplying the predetermined matrix by the second state variable vector by the second state variable vector, and performing probabilistic updating 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 using an inner product of the first state variable vector and the second vector to calculate a function value of the objective function. Effect of the Invention

[0010] 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. [Brief description of the drawings]

[0011] [Figure 1] FIG. 1 is a conceptual diagram illustrating a relationship between a variable array and an objective function value of an optimization problem. [Diagram 2] FIG. 2 is a diagram illustrating an example of the relationship between variables of an objective function of an optimization problem. [Diagram 3] FIG. 2 is a diagram illustrating an example of the relationship between variables of an objective function of an optimization problem. [Figure 4]FIG. 13 is a diagram illustrating an example of an optimal solution search process for multiple replicas. [Diagram 5] FIG. 2 is a diagram illustrating hardware of a computer that realizes the information processing device according to the embodiment. [Figure 6] 1 is a functional block diagram illustrating an arithmetic circuit constituting an arithmetic device according to an embodiment; [Figure 7] 1 is a functional block diagram illustrating a configuration of an information processing apparatus according to an embodiment. [Figure 8] 11 is a flowchart illustrating an optimal solution search process executed by the information processing device according to the embodiment. [Figure 9] 13 is a flowchart illustrating details of an interaction calculation process. [Figure 10] FIG. 11 is a diagram illustrating a simulation result of a transition in the number of replicas due to resampling. [Figure 11] 11A and 11B are diagrams illustrating simulation results comparing the calculation speed according to the embodiment with a comparative example. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0012] Hereinafter, an embodiment of the present invention will be described with reference to the drawings. The embodiment is an example for explaining the present invention, and appropriate omissions and simplifications are made for clarity of explanation. The present invention can be implemented in various other forms. Unless otherwise specified, each component may be singular or plural.

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

[0014] In the embodiment, a process performed by executing a program may be described. Here, a computer executes the program using a processor (e.g., CPU, GPU), and performs the process defined by the program while using a storage resource (e.g., memory) and an interface device (e.g., communication port), etc. Therefore, the subject of the process performed by executing the program may be the processor. Similarly, the subject of the process performed by executing the program may be a controller, device, system, computer, or node having a processor. The subject of the process performed by executing the program may be a calculation unit, and may include a dedicated circuit that performs a specific process. 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), etc.

[0015] The program may be installed in 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 for storing 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.

[0016] (Theoretical Background of the Embodiment) Let s1,…,s be the variables of the optimization problem that optimizes the objective function associated with the interaction model. N Let N be the number of variables, and each variable s i Domain D of i Let be either the binary value {-1, +1} or the continuous value [-1, +1]. i Domain D of iWhich is the case is determined for each problem. The objective function H of the optimization problem is expressed as Equation 1. In other words, the objective function H is expressed as a quadratic function of the variable s.

[0017]

number

[0018]

number

[0019]

number

[0020]

number

[0021] If for all i, D i ={-1, +1}, 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 by an algorithm using the Markov Chain Monte Carlo method (hereinafter referred to as MCMC).

[0022] 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' that is 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.

[0023] 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 the optimal solution or 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)) proposed in Non-Patent Document 1.

[0024] Instead of solving the minimization problem shown in Equation 4, consider solving the minimization problem shown in Equation 5, where S relaxed ={s|s i ∈[-1,+1]}.

[0025]

number

[0026] The optimal solution of Eq. 5 is * =[s1 * ,…,s N * ]. We will not go into the proof here, but the s obtained from Eq. 6 is + =[s1 + ,…,s N + ] is one of the optimal solutions of Equation 4. The goal of this embodiment is to find the optimal solution of Equation 2, but the optimal solution s* Even if we obtain the transformation of Eq. 6 after solving + Here, the function sgn returns +1 if the argument is 0 or greater, and -1 otherwise.

[0027]

number

[0028] Here, N-dimensional vector v=[v 1, …,v N ] is introduced to define the function H' shown in Equation 7.

[0029]

number

[0030] Here, the function V(v) is defined as shown in Equation 8.

[0031]

number

[0032] Matrix W=diag(w1,…,w N ) is an arbitrary diagonal matrix, and v i is a real number that ranges from [-1, +1]. Instead of the minimization problem shown in Equation 5, we introduce Equation 9, which is the minimization problem of H'(s, v).

[0033]

number

[0034] Define two N-dimensional vectors, x = s + v and y = sv. The objective function of the optimization problem to be solved is originally only H, but by introducing a function called V, it becomes possible to obtain a new function that can be updated in parallel using MCMC. Then, the function H' can be rewritten as Equation 10.

[0035]

number

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

[0037]

number

[0038] The optimal solution of Equation 11 is x * , y * Then, s * =(x * +y * ) / 2. These arguments also hold when W is a zero matrix.

[0039] From the above, the optimal solution to the optimization problem expressed by Equation 2 can be obtained from the solution to the constrained quadratic programming problem shown in Equation 11. MCMC is used to obtain this solution.

[0040] FIG. 2 is a diagram illustrating the relationship between each variable of the objective function of an optimization problem, and shows a graphical model showing the relationship between each variable of function G in Equation 11. FIG. 2 illustrates the case of N=6, but the same applies for general N. The relationship between each variable of function G can be expressed as a complete bipartite graph. If there is a variable x i The variables multiplied by are y1,…,y N and x i When MCMC stochastically updates a variable value, it uses the values ​​of variables related to that variable. In other words, when updating the value of variable x1, y1,…,y N and x1, and the other variables (here, x2,…,x N ) is not referenced. This is also true for updating other variables, such as the value of x2. Therefore, if the value of the variable array y is constant, the theoretical requirements of MCMC are not violated even if each value of the array x is independently and simultaneously updated stochastically.

[0041] Similarly, the variable y i The variables multiplied by are also x1,…,xN andy i Therefore, when the value of the variable array x is constant, each value of the array y can be independently and simultaneously updated probabilistically.

[0042] From the above, "x1,...,x N Simultaneous update of y1,…,y N By executing MCMC, which is composed of a procedure of repeating "simultaneous update of G", it is possible to search for arrays x and y that minimize the function G while taking advantage of the speedup achieved by parallelization.

[0043] It should be noted that in the discussion of this embodiment, no constraints are placed on the matrix J. For example, even if all elements of the matrix J are nonzero, the above discussion still holds, and parallel updating is possible.

[0044] Fig. 3 is an example of a fully connected graph, illustrating the relationship between the variables of the objective function of an optimization problem. Fig. 3 illustrates the case of N=6, but the same applies to 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 expressed by a fully connected graph as shown in Fig. 3, so that the probability of only one variable at a time can be updated, and updates are limited to sequential updates.

[0045] From here, we will describe the probabilistic update procedure for each variable. Let x i Let variables y1,…,y N For a given value of , the variable x in the Boltzmann distribution at temperature T i The probability of existence of p(x i ) satisfies Equation 12.

[0046]

number

[0047] However, variable A i is the value calculated using Equation 13.

[0048]

number

[0049] Variable x i andy i is |x i |+|y i |≦2, so x i The range in which it can be moved is -(2-|y i |) or more (2-|y i |) or less. Therefore, the variable x i is average A i / w i , distributed T / w i In the normal distribution of -(2-|y i |) or more (2-|y i Based on the truncated normal distribution with a domain of |), x i In this method, we can sample the next state of x i This means that the next state is determined regardless of the current state of y i In this specification, when there is no distinction between the variables x and y, they may be expressed as s.

[0050] Random numbers following the standard normal distribution can be generated by the Box-Muller method. Since the domain of definition is limited, the algorithm shown in Non-Patent Document 2 can be used.

[0051] The search for an optimal solution can be regarded as sampling from an equilibrium state at temperature 0. Therefore, in order to achieve a high-quality solution search, it is preferable to converge to the equilibrium state in a short time. In order to improve the convergence to the 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. This samples not just one, but K states from the Boltzmann distribution of temperature T as candidates for the next state. Then, in addition to the sampled K states, a total of (K+1) states of the current state are rearranged to obtain x c 0 <=…<=x c r =x i ≦x c KIn other words, the current state is the (r+1)th smallest of the (K+1) values. And x c K+1-r is adopted as the next state. In this method, the next state is x i It depends on the current state of

[0052] As described above, the optimal solution search is performed by a stochastic search process that repeats stochastic transitions of variables x and y that represent the state of the interaction model. However, as a result of the stochastic transition, there are cases where the optimal solution or an approximate solution cannot be obtained even after many trials, as a result of falling into a local solution. 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 update without falling into a local solution and speeding up the optimal solution search.

[0053] 4 is a diagram illustrating an example of an optimal solution search process for multiple replicas, in which one replica is placed in each node.

[0054] For example, multiple nodes with their initial states set at temperature T = T0 (time step 0), for example replicas 1-1, 1-2, and 1-3 placed on three nodes #1, #2, and #3, are used to independently and in parallel search for optimal solutions. t At 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 searched for optimal solutions up until that point, are calculated. The objective function of the j-th replica 1-j (j = 1, , R) among the R replicas is calculated as function G j Let us assume that.

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

[0056]

number

[0057] where Δβ t is the temperature T at time step t when resampling t Inverse temperature β t (=1 / T t ) difference (β t -β t-1 ) The weight of replica 1-j is α j From the above, the occurrence probability q j is calculated as in Equation 15.

[0058]

number

[0059] However, the denominator of the right-hand side of Equation 15 is the weights α of all replicas 1-j before resampling. j That is, Equation 15 expresses the sum of the weight α j Normalize.

[0060] Then, as the replicas after resampling, the number of replicas selected for each replica 1-j is N j is determined as in Equation 16. The selection number N j is the number of replicas 1-j used in the search for the optimal solution at the next time step. j The decimal places of R are rounded appropriately. R represents the total number of replicas, which is constant before and after resampling, but may increase or decrease after resampling.

[0061]

number

[0062] And the number of selections N j The destination of the replica after resampling is determined according to, and the information of each replica 1-j is copied as necessary.

[0063] In the example of FIG. 4, R=3. Replica 1-1 is R×q1=0.6, replica 1-2 is R×q2=2.1, and replica 1-3 is R×q3=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 probabilistic search process is executed in parallel and independently for each of one replica 1-1 and two replicas 1-2.

[0064] In the search for the optimal solution, the replicas are executed in parallel, and when an execution termination condition, such as the temperature dropping to a specified temperature, is satisfied, an optimal solution or an approximate solution is output based on variables representing the state of the replica among the multiple replicas for which the value of function G is the smallest.

[0065] Based on the above, the configuration of the information processing device 10 according to this embodiment will be described with reference to FIGS.

[0066] (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, by using a virtual information processing resource such as a cloud server, a part or whole of which is provided 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.

[0067] The processor 11 is configured using, for example, a CPU (Central Processing Unit) or an MPU (Micro Processing Unit). The main storage 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 storage device 13 is a HDD (Hard Disk Drive), a flash memory, SSD (Solid State Drive), an optical storage device (CD (Compact Disc), DVD (Digital Versatile Disc), etc.), etc. The programs and data stored in the auxiliary storage device 13 are read into the main storage device 12 as needed.

[0068] The input device 14 is a user interface that accepts information input from a 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 a user, and is, for example, a display device (LCD (Liquid Crystal Display), a graphic 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.

[0069] The arithmetic device 17 is a device that executes processing related to an optimal solution search for a combinatorial optimization problem. The arithmetic device 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 device 17 is configured by hardware such as a CMOS (Complementary Metal Oxide Semiconductor) circuit, FPGA, ASIC, etc.

[0070] The arithmetic device 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 device 17 may be communicably connected to another arithmetic device 17 via a communication line, for example, and may operate in cooperation with the other arithmetic device 17. The functions realized by the arithmetic device 17 may be realized, for example, by having a processor (CPU, GPU, etc.) execute a program.

[0071] 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 acquired 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.

[0072] The information processing device 10 may be an information processing system in which a plurality of devices cooperate to execute processing. A program for realizing the information processing system causes each device to realize each function by each program, thereby realizing the same function as the information processing device 10.

[0073] (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 computer, etc., for example, a difference calculation block 514, a sampling block 515, a next state determination block 516, etc., which will be described later. However, the present invention is not limited to this, and may be implemented by hardware. The arithmetic circuit 500 is implemented by software using a number of arithmetic units mounted on a GPU, vector computer, etc., for example, a difference calculation block 514, a sampling block 515, a next state determination block 516, etc., which will be described later. However, the present invention is not limited to this, and may be implemented by hardware. N or variable array y1,…,y N from the Boltzmann distribution (Equation 12) at temperature T.

[0074] 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 .

[0075] The variable memory 511 of each arithmetic circuit 500 stores the above-mentioned variables x1, . . . , x N and y1,…,y N Information indicating this is stored (see Figure 2).

[0076] The nonlinear coefficient memory 512 stores information representing a 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 a vector h.

[0077] The arithmetic unit 17 receives a control signal EN, a weight signal SW, and a temperature signal TE.

[0078] 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 controls the variables x1, ..., x Nand simultaneously update the variables y1,…,y N 5, for the sake of simplicity, the control signal EN is input only to the sampling block 515, but it is also input to other parts such as the variable memory 511 which require this signal.

[0079] The weighting signal SW is a signal representing a vector of N elements representing the diagonal elements of the diagonal matrix W.

[0080] The difference calculation block 514 receives the values ​​of the matrix J stored in the nonlinear coefficient memory 512, the vector h stored in the linear coefficient memory 513, the weight signal SW, and the variable s (x or y) stored in the variable memory 511. When the control signal EN is H, the difference calculation block 514 calculates (J+diag(w1, . . . , w N ))y+h, when EN is L, (J+diag(w1, ,w N ))x+h. This output value is the A i is equivalent to.

[0081] The sampling block 515 receives the output of the difference calculation block 514, the weight signal SW, the temperature signal TE that holds the value of the temperature parameter, the control signal EN, and the values ​​of other variables. Then, as the i-th element, when the control signal EN is H, −(2−|y i |) or more (2-|y i |) below, when EN is L, -(2-|x i |) or more (2-|x i The output is a random sample from the truncated normal distribution expressed by Equation 12, with the domain being less than or equal to |.

[0082] 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 a simple heat bath method is set as the update rule of the MCMC, the next state determination block 516 receives only one output value of the sampling block 515 and writes it directly into the variable memory 511. If the well-known over-relaxation method is used as the update rule of the MCMC, 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 into the variable memory 511. As is well known, in the over-relaxation method, the next state is determined so that the correlation with the immediately preceding state is negative.

[0083] (Functional configuration of information processing device 10) 7 is a functional block diagram illustrating the configuration of an information processing device 10 according to an 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.

[0084] 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 reading and executing a program stored in the auxiliary storage unit 13 in cooperation with the main storage unit 12, or by hardware included in the calculation unit 17. In addition to the above functions, the information processing device 10 may have other functions such as an operating system, a file system, a device driver, and a DBMS (Data-Base Management System).

[0085] The storage unit 102 is realized by the main storage device 12 or the auxiliary storage 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.).

[0086] The quadratic programming format 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 into the domain data 102c. The domain data 102c indicates, for example, whether each variable takes a binary value or a real 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.

[0087] 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 need to be implemented in the information processing device 10, and the information processing device 10 may import the quadratic programming problem data 102b generated by another information processing device or the like via the input device 14 or the communication device 16.

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

[0089] 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.

[0090] The variable value initialization unit 101d initializes the value of each variable stored in the variable memory 511 of the arithmetic circuit 500. For example, the variable value initialization unit 101d may randomly sample and determine the value of each variable uniformly from -1 to +1. At this time, the constraint on the variables, |x i |+|y i Care must be taken to ensure that |≦2 is satisfied. Also, please note that the values ​​of each variable at this time are treated as continuous values.

[0091] The temperature setting unit 101e sets the temperature T used when the interaction calculation executing unit 101f searches for an optimal solution.

[0092] 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 a function G expressed by Equation 11 for each temperature T set by the temperature setting unit 101e. In the interaction calculation process, the interaction calculation execution unit 101f changes the temperature T, for example, from a higher temperature to a lower temperature.

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

[0094] In addition, the interaction calculation execution unit 101f executes resampling at a predetermined period. That is, the interaction calculation execution unit 101f calculates the energy (objective function value (function Gj Calculate the value of

[0095] Here, the function G j The multiplication of the matrix J and the variable y (the first term on the right side of Equation 11) is a heavy calculation load. However, this multiplication of the matrix J and the variable y was performed in the previous update step of the variable x, and the multiplication result Jy is stored in the work area. Therefore, the function G j When calculating the value of J, the multiplication result Jy stored in the work memory is used instead of multiplying the matrix J by the variable y. Since the multiplication result Jy is a vector, the function G j The first term on the right hand side of the formula is 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 amount of calculation of the objective function value can be suppressed, and the calculation speed of the optimization search can be increased.

[0096] Then, the interaction calculation execution unit 101f calculates the number of selected replicas 1-j after resampling, N, from Equations 14, 15, and 16 based on the energy of each replica 1-j at time step t. j Determine.

[0097] Then, the interaction calculation execution unit 101f selects the number N of each replica 1-j. j Among them, the copy destination of each replica 1-j of the increased portion after resampling is determined. 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 that 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 the respective upper limits. 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.

[0098] In this way, by giving priority to copying 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.

[0099] 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 to equation 11. In accordance with the above discussion, it calculates the N-dimensional vector s*=(x+y) / 2. Then it reads out the domain data 102c, and outputs the vector s+ obtained by equation 6 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}, it outputs sgn(s*i), and if the i-th domain is [-1, +1], it outputs s i In this way, a solution is found according to the defined range.

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

[0101] 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 by 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. In addition, 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).

[0102] The processes of S12 to S15 are executed for R replicas created by R arithmetic circuits 500 of the arithmetic device 17.

[0103] 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 by, for example, the input device 14, the output device 15, the communication device 16, etc.).

[0104] In S13, the weight setting unit 101c determines the value of the weight signal SW. As explained in the above-mentioned formula 8, the weight signal SW is allowed to take any value in 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 formulas 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 executed in the arithmetic unit 17 or the processor 11. Alternatively, the user may set it by himself (S13).

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

[0106] In S15, the temperature setting unit 101e sets a series T γ (γ=1, 2, 3, . . . ) is initially set. Note that the above subscript γ indicates the type of temperature T to be set. The method of setting the temperature T may be, for example, the method disclosed in Patent Document 1.

[0107] In S16, the interaction calculation execution unit 101f executes the interaction calculation process, the details of which will be described later with reference to FIG.

[0108] In S17, the interaction calculation executing unit 101f judges whether or not the stop condition is satisfied (for example, whether or not the temperature T has reached a preset minimum temperature). If the interaction calculation executing unit 101f judges that the stop condition is satisfied (S17: Yes), the process proceeds to S18, whereas if the interaction calculation executing unit 101f judges that the stop condition is not satisfied (S17: No), the process returns to S16.

[0109] In S18, the variable value reading unit 101g reads the function G j The variable value reading unit 101g reads out the variable value stored in the variable memory 511 from the arithmetic circuit 500 of replica 1-j having the smallest value of . The variable value reading unit 101g also reads out 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.

[0110] (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 the calculations of each calculation circuit 500.

[0111] In S162, the interaction calculation execution unit 101f calculates the function G j Calculate the value of (Equation 11). As described above, the function G j When calculating the objective function value of the function G, the calculation results of the update process of the variables x and y are used. j The calculation of the value of may be performed by the arithmetic circuit 500 of each replica 1-j.

[0112] In S163, the interaction calculation execution unit 101f selects the number N j And determine the destination to which the incremental replica will be copied.

[0113] In S164, the interaction calculation executing unit 101f transmits information about the replica to be copied to the copy destination node (the own node or another node) determined in S163.

[0114] (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 that have the same original replica at time 0 are all in the same lineage. 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.

[0115] (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 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.

[0116] Furthermore, in the disclosed technology, state updating and function value calculation are performed independently and in parallel for multiple replicas corresponding to the first state variable vector and the second state variable vector, so that state updating can be performed globally without being limited to individual local solutions.

[0117] Furthermore, in the disclosed technology, the state update and function value calculation for each replica are executed by a plurality of arithmetic units, so that parallel processing can be performed while distributing the load, thereby increasing the processing speed.

[0118] 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 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.

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

[0120] Furthermore, in the disclosed technology, the number of replicas selected in resampling is greater than the total number of replicas before resampling, so that the optimum solution search can be performed using abundant resources with many system replicas, further speeding up the optimum solution search.

[0121] Furthermore, in the disclosed technology, the number of replicas selected in resampling is smaller than the total number of replicas before resampling, and therefore the search for an optimal solution can be efficiently accelerated by allocating limited resources to a series of replicas with high objective function values.

[0122] Furthermore, in the disclosed technology, if a copy of a replica can be placed on a computing device on which a selected replica is placed according to the number of replicas selected in resampling, the computing device is determined as the copy destination on which the copy of the replica is preferentially placed. On the other hand, if placement is not possible, another computing device is determined as the copy destination. Therefore, the occurrence of inter-node communication is suppressed when transferring (copying) replica information during resampling, which suppresses communication time and speeds up the search for an optimal solution. EXAMPLES

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

[0124] Fig. 11 is a diagram illustrating a simulation result in which the calculation speed according to the embodiment is compared with a comparative example. In Fig. 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 performed for each run with different random number seeds, and the average value, first quartile, and third quartile of the objective function value at each time were plotted. The dashed line indicates the average value, the outer line at the lower end of the gray line overlapping the dashed line indicates the first quartile, and the outer line at the upper end indicates the third quartile.

[0125] 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 GPUs. 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 GPUs.

[0126] [Table 1]

[0127] As can be seen from FIG. 11, Examples 1 and 2 are in the lowest energy state (objective function value=−4.40×10 5 ) faster. Moreover, Example 1 reached near the lowest energy state faster than Comparative Example 2, which also had one GPU implementation. Moreover, Example 2 reached near the lowest energy state faster than Comparative Example 3, which also had four GPU implementation. Moreover, Example 2 reached near the lowest energy state faster than Example 1.

[0128] From the above simulation results, it was found 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.

[0129] (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) and a quantum neural network 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.

[0130] As an implementation example 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 using SRAM (Static Random Access Memory) technology may be arranged, and a memory for storing variables and a circuit for updating the variables may be arranged in each unit.

[0131] The present invention is not limited to the above-described embodiments, and includes various modified examples. 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. In addition, it is also possible to replace a 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, as long as there is no contradiction. In addition, it is possible to add, delete, replace, integrate, or distribute a part of the configuration of each embodiment. In addition, the configurations and processes shown in the embodiments can be appropriately distributed, integrated, or replaced based on processing efficiency or implementation efficiency. [Explanation of symbols]

[0132] 10: information processing device, 11: processor, 12: main memory device, 13: auxiliary memory 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

1. 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 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; An information processing system comprising:

2. 2. The information processing system according to claim 1, The processing unit includes: The state update process and the function value calculation process are 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. An information processing system comprising:

3. 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. An information processing system comprising:

4. 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. An information processing system comprising:

5. 5. The information processing system according to claim 4, The selected number is equal to the total number of the plurality of replicas. An information processing system comprising:

6. 5. The information processing system according to claim 4, The selected number is greater than the total number of the plurality of replicas. An information processing system comprising:

7. 5. The information processing system according to claim 4, The selected number is less than the total number of the plurality of replicas. An information processing system comprising:

8. 5. The information processing system according to claim 4, The processing unit includes: In the resampling process, if a copy of the replica can be placed in the arithmetic device on which the selected replica is placed, the arithmetic device is determined as a copy destination for placing the copy of the replica, depending on the selection number, and if the copy cannot be placed in the arithmetic device, another arithmetic device is determined as the copy destination. An information processing system comprising:

9. 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; 2. An information processing method comprising:

10. 10. The information processing method according to claim 9, The state update process and the function value calculation process include: 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.

23. An information processing method comprising:

11. 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.

23. An information processing method comprising:

12. 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; 2. An information processing method comprising:

13. 13. The information processing method according to claim 12, In the resampling process, if the arithmetic device on which the selected replica is to be placed can place a copy of the replica, the arithmetic device is determined as a copy destination for placing the copy of the replica, depending on the selection number, and if the arithmetic device cannot place the copy, another arithmetic device is determined as the copy destination.

23. An information processing method comprising:

14. A program for causing a computer to function as the information processing system according to any one of claims 1 to 8.

Images