WO2020196872A1 - 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム - Google Patents

情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム Download PDF

Info

Publication number
WO2020196872A1
WO2020196872A1 PCT/JP2020/014190 JP2020014190W WO2020196872A1 WO 2020196872 A1 WO2020196872 A1 WO 2020196872A1 JP 2020014190 W JP2020014190 W JP 2020014190W WO 2020196872 A1 WO2020196872 A1 WO 2020196872A1
Authority
WO
WIPO (PCT)
Prior art keywords
variable
vector
coefficient
exchange
information processing
Prior art date
Application number
PCT/JP2020/014190
Other languages
English (en)
French (fr)
Inventor
鈴木　賢
隼人 後藤
光介 辰村
Original Assignee
株式会社 東芝
東芝デジタルソリューションズ株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社 東芝, 東芝デジタルソリューションズ株式会社 filed Critical 株式会社 東芝
Priority to CA3135142A priority Critical patent/CA3135142A1/en
Priority to CN202080025407.3A priority patent/CN113646787A/zh
Priority to JP2021509666A priority patent/JP7421545B2/ja
Publication of WO2020196872A1 publication Critical patent/WO2020196872A1/ja
Priority to US17/487,179 priority patent/US20220012017A1/en

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00—Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/38—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation
    • G06F7/48—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation using non-contact-making devices, e.g. tube, solid state device; using unspecified devices
    • G06F7/57—Arithmetic logic units [ALU], i.e. arrangements or devices for performing two or more of the operations covered by groups G06F7/483 – G06F7/556 or for performing logical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N10/00—Quantum computing, i.e. information processing based on quantum-mechanical phenomena
    • G06N10/60—Quantum algorithms, e.g. based on quantum optimisation, quantum Fourier or Hadamard transforms

Definitions

  • An embodiment of the present invention relates to an information processing device, an information processing system, an information processing method, a storage medium, and a program.
  • the combinatorial optimization problem is a problem of selecting the most suitable combination from a plurality of combinations.
  • Combinatorial optimization problems are mathematically reduced to a problem called an "objective function" that maximizes a function having a plurality of discrete variables or minimizes the function.
  • the combinatorial optimization problem is a universal problem in various fields such as finance, logistics, transportation, design, manufacturing, and life science, but the number of combinations increases in the order of the exponential function of the problem size, so-called "combinatorial explosion". Therefore, it is not always possible to find the optimum solution. Moreover, it is often difficult to even obtain an approximate solution close to the optimum solution.
  • An embodiment of the present invention provides an information processing device, an information processing system, an information processing method, a storage medium, and a program for calculating a solution of a combinatorial optimization problem within a practical time.
  • the information processing apparatus is configured to repeatedly update the first vector having the first variable as an element and the second vector having the second variable corresponding to the first variable as an element. It is equipped with a plurality of arithmetic circuits and a data exchange circuit. Each of the arithmetic circuits updates the first variable based on the corresponding second variable, weights the first variable with a first coefficient, adds it to the corresponding second variable, and a plurality of the first variables. The problem term is calculated using the above, and the problem term is added to the second variable. In each of the arithmetic circuits, different values are set for the first coefficient, and the data exchange circuit is used. Is configured to perform at least one of the exchange of the first vector and the second vector, or the exchange of the first coefficient between the arithmetic circuits.
  • the figure which showed the configuration example of an information processing system A block diagram showing a configuration example of the management server.
  • the figure which shows the example of the data stored in the storage of the calculation server. A flowchart showing an example of processing when calculating the solution of a simulated branch algorithm by time evolution.
  • a flowchart showing an example of processing when a calculation is performed using a fluctuating coefficient p (m) A flowchart showing an example of processing in the case of performing replica exchange in parallel.
  • FIG. 1 is a block diagram showing a configuration example of the information processing system 100.
  • the information processing system 100 of FIG. 1 includes a management server 1, a network 2, calculation servers (information processing devices) 3a to 3c, cables 4a to 4c, a switch 5, and a storage device 7. Further, FIG. 1 shows a client terminal 6 capable of communicating with the information processing system 100.
  • the management server 1, the calculation servers 3a to 3c, the client terminal 6, and the storage device 7 can communicate data with each other via the network 2.
  • the calculation servers 3a to 3c can store data in the storage device 7 and read data from the storage device 7.
  • the network 2 is, for example, the Internet in which a plurality of computer networks are connected to each other.
  • the network 2 can use wired, wireless, or a combination thereof as a communication medium.
  • an example of the communication protocol used in the network 2 is TCP / IP, but the type of the communication protocol is not particularly limited.
  • the calculation servers 3a to 3c are connected to the switch 5 via cables 4a to 4c, respectively. Cables 4a-4c and switch 5 form an interconnect between compute servers. The calculation servers 3a to 3c can also perform data communication with each other via the interconnect.
  • the switch 5 is, for example, an Infiniband switch.
  • the cables 4a to 4c are, for example, Infiniband cables. However, a wired LAN switch / cable may be used instead of the InfiniBand switch / cable.
  • the communication standard and communication protocol used in the cables 4a to 4c and the switch 5 are not particularly limited. Examples of the client terminal 6 include a notebook PC, a desktop PC, a smartphone, a tablet, an in-vehicle terminal, and the like.
  • the management server 1 controls the calculation server by converting, for example, the combinatorial optimization problem input by the user into a format that can be processed by each calculation server. Then, the management server 1 acquires the calculation result from each calculation server and converts the aggregated calculation result into a solution of the combinatorial optimization problem. In this way, the user can obtain a solution to the combinatorial optimization problem.
  • the solution of the combinatorial optimization problem shall include an optimal solution and an approximate solution close to the optimal solution.
  • FIG. 1 shows three calculation servers.
  • the number of calculation servers included in the information processing system is not limited.
  • the number of calculation servers used to solve the combinatorial optimization problem is not particularly limited.
  • the information processing system may include one calculation server.
  • the combinatorial optimization problem may be solved by using one of a plurality of calculation servers included in the information processing system.
  • the information processing system may include hundreds or more calculation servers.
  • the calculation server may be a server installed in a data center or a desktop PC installed in an office. Further, the calculation server may be a plurality of types of computers installed at different locations.
  • the type of information processing device used as the calculation server is not particularly limited.
  • the calculation server may be a general-purpose computer, a dedicated electronic circuit, or a combination thereof.
  • FIG. 2 is a block diagram showing a configuration example of the management server 1.
  • the management server 1 in FIG. 2 is, for example, a computer including a central processing unit (CPU) and a memory.
  • the management server 1 includes a processor 10, a storage unit 14, a communication circuit 15, an input circuit 16, and an output circuit 17. It is assumed that the processor 10, the storage unit 14, the communication circuit 15, the input circuit 16 and the output circuit 17 are connected to each other via the bus 20.
  • the processor 10 includes a management unit 11, a conversion unit 12, and a control unit 13 as internal components.
  • the processor 10 is an electronic circuit that executes calculations and controls the management server 1.
  • the processor 10 for example, a CPU, a microprocessor, an ASIC, an FPGA, a PLD, or a combination thereof can be used.
  • the management unit 11 provides an interface for operating the management server 1 via the user's client terminal 6. Examples of the interface provided by the management unit 11 include API, CLI, and a web page.
  • the user can input information on the combinatorial optimization problem via the management unit 11, and can view and / or download the calculated solution of the combinatorial optimization problem.
  • the conversion unit 12 converts the combinatorial optimization problem into a format that can be processed by each calculation server.
  • the control unit 13 transmits a control command to each calculation server.
  • control unit 13 After the control unit 13 acquires the calculation result from each calculation server, the conversion unit 12 aggregates the plurality of calculation results and converts them into a solution of the combinatorial optimization problem. Further, the control unit 13 may specify the processing content to be executed by each calculation server or the processor in each server.
  • the storage unit 14 stores various data including the program of the management server 1, the data necessary for executing the program, and the data generated by the program.
  • the program includes both the OS and the application.
  • the storage unit 14 may be a volatile memory, a non-volatile memory, or a combination thereof. Examples of volatile memory include DRAM, SRAM and the like. Examples of non-volatile memory include NAND flash memory, NOR flash memory, ReRAM, or MRAM. Further, as the storage unit 14, a hard disk, an optical disk, a magnetic tape, or an external storage device may be used.
  • the communication circuit 15 transmits / receives data to / from each device connected to the network 2.
  • the communication circuit 15 is, for example, a wired LAN NIC (Network Interface Card). However, the communication circuit 15 may be another type of communication circuit such as a wireless LAN.
  • the input circuit 16 realizes data input to the management server 1. It is assumed that the input circuit 16 includes, for example, USB, PCI-Express, or the like as an external port.
  • the operating device 18 is connected to the input circuit 16.
  • the operation device 18 is a device for inputting information to the management server 1.
  • the operating device 18 is, for example, a keyboard, a mouse, a touch panel, a voice recognition device, and the like, but is not limited thereto.
  • the output circuit 17 realizes data output from the management server 1. It is assumed that the output circuit 17 is provided with HDMI, DisplayPort, or the like as an external port.
  • the display device 19 is connected to the output circuit 17. Examples of the display device 19 include, but are not limited to, an LCD (liquid crystal display), an organic EL (organic electroluminescence) display, or a projector.
  • the administrator of the management server 1 can perform maintenance on the management server 1 by using the operation device 18 and the display device 19.
  • the operation device 18 and the display device 19 may be incorporated in the management server 1. Further, the operation device 18 and the display device 19 do not necessarily have to be connected to the management server 1. For example, the administrator may perform maintenance on the management server 1 using a client terminal capable of communicating with the network 2.
  • FIG. 3 shows an example of data stored in the storage unit 14 of the management server 1.
  • the problem data 14A, the calculation data 14B, the management program 14C, the conversion program 14D, and the control program 14E are stored in the storage unit 14 of FIG.
  • the problem data 14A includes data of a combinatorial optimization problem.
  • the calculation data 14B includes the calculation results collected from each calculation server.
  • the management program 14C is a program that realizes the functions of the management unit 11 described above.
  • the conversion program 14D is a program that realizes the functions of the conversion unit 12 described above.
  • the control program 14E is a program that realizes the functions of the control unit 13 described above.
  • FIG. 4 is a block showing a configuration example of the calculation server.
  • FIG. 4 illustrates the configuration of the calculation server 3a as an example.
  • the other calculation server may have the same configuration as the calculation server 3a, or may have a configuration different from that of the calculation server 3a.
  • the calculation server 3a is, for example, an information processing device that executes the calculation of the first vector and the second vector independently or in sharing with other calculation servers.
  • at least one of the calculation servers may calculate the problem terms between the elements of the first vector.
  • the problem term may be calculated based on the Ising model described later.
  • the problem section may include many-body interactions described later.
  • m indicates an integer replica number from 1 to M. The details of the replica will be described later.
  • the calculation server 3a includes, for example, a communication circuit 31, a shared memory 32, processors 33A to 33D, a storage 34, and a host bus adapter 35. It is assumed that the communication circuit 31, the shared memory 32, the processors 33A to 33D, the storage 34, and the host bus adapter 35 are connected to each other via the bus 36.
  • the communication circuit 31 transmits / receives data to / from each device connected to the network 2.
  • the communication circuit 31 is, for example, a wired LAN NIC (Network Interface Card). However, the communication circuit 31 may be another type of communication circuit such as a wireless LAN.
  • the shared memory 32 is a memory that can be accessed from the processors 33A to 33D. Examples of the shared memory 32 include volatile memories such as DRAM and SRAM. However, as the shared memory 32, another type of memory such as a non-volatile memory may be used.
  • the processors 33A to 33D can share data via the shared memory 32. It should be noted that not all the memories of the calculation server 3a need to be configured as shared memory. For example, a part of the memory of the calculation server 3a may be configured as a local memory that can be accessed only by any processor.
  • Processors 33A to 33D are electronic circuits that execute calculation processing.
  • the processor may be, for example, a CPU (Central Processing Unit), a GPU (Graphics Processing Unit), an FPGA (Field-Programmable Gate Array), or an ASIC (Application Specific Integrated Circuit), which may be a combination of these or a combination of these. You may.
  • the processor may be a CPU core or a CPU thread.
  • the processor may be connected to other components of the calculation server 3a via a bus such as PCI express.
  • the calculation server is equipped with four processors.
  • the number of processors provided in one computing server may be different from this.
  • the number and / or type of processors implemented in the compute server may be different.
  • the storage unit of the information processing apparatus may be configured to store the first variable which is an element of the first vector and the second variable which is an element of the second vector.
  • the plurality of arithmetic circuits 50 are configured to repeatedly update the first vector whose element is the first variable and the second vector whose element is the second variable corresponding to the first variable. For example, each arithmetic circuit 50 updates the first variable based on the corresponding second variable, weights the first variable with the first coefficient, adds it to the corresponding second variable, and uses a plurality of first variables. It is configured to calculate the problem term and add the problem term to the second variable. Further, different values may be set for the first coefficient in each arithmetic circuit 50.
  • the data exchange circuit 51 is configured to perform at least one of the exchange of the first vector and the second vector, or the exchange of the first coefficient between the arithmetic circuits 50.
  • the problem terms calculated by the arithmetic circuit may be based on the Ising model. In addition, the problem terms calculated by the arithmetic circuit may include many-body interactions. The details of the problem items will be described later.
  • the processors 33A to 33C correspond to a plurality of arithmetic circuits
  • the processor 33D corresponds to a data exchange circuit.
  • the correspondence between the arithmetic circuit / data exchange circuit shown in FIG. 4 and the processor is only an example. Therefore, the correspondence between the arithmetic circuit / data exchange circuit and the processor may be different from this.
  • the number of processors assigned to the arithmetic circuit / data exchange circuit is not particularly limited. As will be described later, the same processor may combine the roles of an arithmetic circuit and a data exchange circuit.
  • a plurality of types of processors for example, CPU, GPU, FPGA
  • different types of processors may be assigned to the arithmetic circuit and the data exchange circuit.
  • the processing content is assigned to each processor.
  • the unit of computational resources to which the processing contents are assigned is not limited.
  • the processing contents may be assigned for each computer.
  • the information processing system may include a plurality of first computers and a second computer.
  • the first computer is configured to repeatedly update the first vector having the first variable as an element and the second vector having the second variable corresponding to the first variable as an element.
  • the second computer is configured to perform data exchange between the first computers.
  • the first computer and the second computer correspond to, for example, the above-mentioned calculation server.
  • Each first computer updates the first variable based on the corresponding second variable, weights the first variable with the first coefficient, adds it to the corresponding second variable, and uses the plurality of first variables to solve the problem. It may be configured to calculate the term and add the problem term to the second variable.
  • the second computer may be configured to perform at least one of the exchange of the first vector and the second vector, or the exchange of the first coefficient between the first computers. Further, different values may be set for the first coefficient in each first computer.
  • the storage 34 stores various data including the program of the calculation server 3a, the data necessary for executing the program, and the data generated by the program.
  • the program includes both the OS and the application.
  • the storage 34 may be a volatile memory, a non-volatile memory, or a combination thereof. Examples of volatile memory include DRAM, SRAM and the like. Examples of non-volatile memory include NAND flash memory, NOR flash memory, ReRAM, or MRAM. Further, as the storage 34, a hard disk, an optical disk, a magnetic tape, or an external storage device may be used.
  • the host bus adapter 35 realizes data communication between calculation servers.
  • the host bus adapter 35 is connected to the switch 5 via a cable 4a.
  • the host bus adapter 35 is, for example, an HCA (Host Channel Adapter).
  • HCA Hyper Channel Adapter
  • FIG. 5 shows an example of data stored in the storage of the calculation server.
  • the calculation data 34A, the calculation program 34B, and the control program 34C are stored in the storage 34 of FIG.
  • the calculation data 34A includes data in the middle of calculation or a calculation result of the calculation server 3a. At least a part of the calculated data 34A may be stored in a different storage hierarchy such as the shared memory 32, the cache of the processor, or the register of the processor.
  • the calculation program 34B is a program that realizes a calculation process in each processor and a data storage process in the shared memory 32 and the storage 34 based on a predetermined algorithm.
  • the control program 34C is a program that controls the calculation server 3a based on the command transmitted from the control unit 13 of the management server 1 and transmits the calculation result of the calculation server 3a to the management server 1.
  • An Ising machine is an example of an information processing device used to solve a combinatorial optimization problem.
  • the Ising machine is an information processing device that calculates the energy of the ground state of the Ising model.
  • the Ising model has often been used mainly as a model for ferromagnets and phase transition phenomena.
  • the Ising model has been increasingly used as a model for solving combinatorial optimization problems.
  • the following equation (1) shows the energy of the Ising model.
  • s i and s j are spins
  • spin is a binary variable having a value of either +1 or -1.
  • N is the number of spins.
  • h i is the local magnetic field acting on each spin.
  • J is a matrix of coupling coefficients between spins.
  • the matrix J is a real symmetric matrix having a diagonal component of 0. Therefore, J ij indicates the elements of the matrix J in row i and column j.
  • the Ising model of Eq. (1) is a quadratic equation for spin, but an extended Ising model (Ising model with many-body interaction) including terms of the third order or higher of spin is used. May be good.
  • the solution of the Ising model is expressed in the form of spin vectors (s 1 , s 2 , ..., S N ).
  • the vector (s 1 , s 2 , ..., S N ) at which the energy E Ising is the minimum value is called an optimum solution.
  • the calculated Ising model solution does not necessarily have to be an exact optimal solution.
  • an Ising problem the problem of finding an approximate solution in which the energy E Ising is as small as possible (that is, an approximate solution in which the value of the objective function is as close to the optimum value as possible) using the Ising model is referred to as an Ising problem.
  • the spin s i of the equation (1) is a binary variable, it can be easily converted with the discrete variable (bit) used in the combinatorial optimization problem by using the equation (1 + s i ) / 2. .. Therefore, it is possible to find the solution of the combinatorial optimization problem by converting the combinatorial optimization problem into the Ising problem and letting the Ising machine perform the calculation.
  • the problem of finding a solution that minimizes the quadratic objective function with a discrete variable (bit) that takes either 0 or 1 as a variable is QUABO (Quadratic Unconstrained Binary Optimization, unconstrained binary variable quadratic optimization). ) It is called a problem. It can be said that the Ising problem represented by the equation (1) is equivalent to the QUABO problem.
  • Quantum annealing realizes quantum annealing using a superconducting circuit.
  • the coherent Ising machine utilizes the oscillation phenomenon of a network formed by an optical parametric oscillator.
  • the quantum branching machine utilizes the quantum mechanical branching phenomenon in a network of parametric oscillators with the Kerr effect. While these hardware implementations have the potential to significantly reduce calculation time, they also have the problem of being difficult to scale up and operate stably.
  • SA Simulated annealing
  • Technology is being developed to perform simulated annealing faster.
  • general simulated annealing is a sequential update algorithm in which each variable is sequentially updated, it is difficult to speed up the calculation process by parallelization.
  • a simulated branching algorithm has been proposed that can solve large-scale combinatorial optimization problems at high speed by parallel calculation in a digital computer.
  • information processing devices, information processing systems, information processing methods, storage media, and programs that solve combinatorial optimization problems using a simulated branching algorithm will be described.
  • H is a Hamiltonian of the following formula (3).
  • each coefficient will be described on the assumption that the time evolution is calculated.
  • the simulated branching algorithm does not necessarily have to calculate the time evolution.
  • the simulated branching algorithm may be solved by the replica exchange method (replica exchange method / parallel tempering) as described later.
  • the coefficient D corresponds to detuning.
  • the coefficient p (t) corresponds to the pumping amplitude.
  • the initial value of the coefficient p (t) may be set to 0.
  • the coefficient K corresponds to a positive Kerr coefficient.
  • a constant coefficient can be used as the coefficient c.
  • the value of the coefficient c may be determined before performing the calculation by the simulated branching algorithm.
  • the coefficient c can be set to a value close to the reciprocal of the maximum eigenvalue of the J (2) matrix.
  • the value c 0.5D ⁇ (N / 2n) can be used.
  • n is the number of edges of the graph related to the combinatorial optimization problem.
  • a (t) is a coefficient that increases with p (t) when calculating the time evolution.
  • ⁇ (p (t) / K) can be used as a (t).
  • the symplectic Euler method can be used to convert (2) above into a discrete recurrence formula to perform the solution.
  • (4) shows an example of a simulated branching algorithm after conversion to a recurrence formula.
  • t is a time
  • ⁇ t is a time step (time step width).
  • the time t and the time step ⁇ t are used to show the correspondence with the differential equation.
  • the time t and the time step ⁇ t do not necessarily have to be included as explicit parameters. For example, if the time step ⁇ t is set to 1, it is possible to remove the time step ⁇ t from the algorithm at the time of mounting.
  • x i (t + ⁇ t) may be interpreted as the updated value of x i (t). That is, "t” in (4) above indicates the value of the variable before the update, and "t + ⁇ t” indicates the value of the variable after the update.
  • the value of spin s i is calculated based on the sign of the variable x i after increasing the value of p (t) from the initial value (eg 0) to a predetermined value.
  • the value of spin s i can be obtained by converting the variable x i with a sign function.
  • Solution of a combinatorial optimization problem e.g., spin s i Ising model
  • the solution (solution vector) of the combinatorial optimization problem may be obtained when the number of updates of the first vector and the second vector or the value of the first coefficient p becomes larger than the threshold value.
  • FIG. 6 shows an example of processing in the case of calculating the solution of the simulated branch algorithm by time evolution. Hereinafter, the process will be described with reference to FIG.
  • the calculation server obtains the matrix J ij and vector h i corresponding to the problem from the management server 1 (step S101). Then, the calculation server initializes the coefficients p (t) and a (t) (step S102). For example, the values of the coefficients p and a can be set to 0 in step S102, but the initial values of the coefficients p and a are not limited.
  • the calculation server initializes the first variable x i and the second variable y i (step S103).
  • the first variable x i is an element of the first vector.
  • the second variable y i is an element of the second vector.
  • the compute server may initialize x i and y i with pseudo-random numbers, for example.
  • the method of initializing x i and y i is not limited. Further, the variables may be initialized at different timings, or at least one of the variables may be initialized a plurality of times.
  • the calculation server updates the first vector by weighting and adding the element y i of the second vector corresponding to the element x i of the first vector (step S104). For example, in step S104, ⁇ t ⁇ D ⁇ y i can be added to the variable x i . Then, the calculation server updates the element y i of the second vector (steps S105 and S106). For example, in step S105, ⁇ t ⁇ [(pdK ⁇ x i ⁇ x i ) ⁇ x i ] can be added to the variable y i . In step S106, it is possible to add a - ⁇ t ⁇ c ⁇ h i ⁇ a - ⁇ t ⁇ c ⁇ ⁇ J ij ⁇ x j and into the variable yi.
  • the calculation server updates the values of the coefficients p and a (step S107). For example, a constant value ( ⁇ p) can be added to the coefficient p, and the coefficient a can be set to the positive square root of the updated coefficient p. However, as will be described later, this is only an example of how to update the values of the coefficients p and a.
  • the calculation server determines whether or not the number of updates of the first vector and the second vector is less than the threshold value (step S108). If the number of updates is less than the threshold value (YES in step S108), the calculation server re-executes the processes of steps S104 to S107.
  • step S109 When the number of updates is equal to or greater than the threshold value (NO in step S108), the spin s i , which is an element of the solution vector, is obtained based on the element x i of the first vector (step S109).
  • step S109 for example, in the first vector, the variable x i is a positive value +1, converts each variable x i is a negative value -1, it is possible to obtain a solution vector.
  • step S108 when the number of updates is less than the threshold value (YES in step S108), the Hamiltonian value may be calculated based on the first vector, and the first vector and the Hamiltonian value may be stored. .. This allows the user to select the approximate solution closest to the optimal solution from the plurality of first vectors.
  • processing may be parallelized using a plurality of calculation servers.
  • Processing may be parallelized by a plurality of processors.
  • the implementation and the mode of parallelization of processing for realizing parallelization of processing are not limited.
  • the execution order of the update processing of the variables x i and y i shown in steps S105 to S106 described above is only an example. Therefore, the update processing of the variables x i and y i may be executed in a different order. For example, the order in which the update process of the variable x i and the update process of the variable y i may be executed may be interchanged. Further, the order of sub-processing included in the update processing of each variable is not limited. For example, the execution order of the addition processing included in the update processing of the variable y i may be different from the example of FIG. The execution order and timing of the processes that are the prerequisites for executing the update process of each variable are not particularly limited.
  • the computing problems terms may be performed in parallel with other processing including processing of updating the variable x i.
  • the order and timing at which the update processing of the variables x i and y i, the sub-processing included in the update processing of each variable, and the calculation processing of the problem item are executed are not limited, and the same applies to the processing of each flowchart shown below. Is.
  • the graph of FIG. 7 shows an example of an update pattern of the coefficient p when calculating the time evolution.
  • the vertical axis of FIG. 7 shows the value of the coefficient p.
  • the horizontal axis of FIG. 7 corresponds to the time or the number of updates.
  • Figure 7 shows that there are multiple update patterns leading to p M is the value of the updated of M times from the initial value p 1. If satisfy the condition at the p 1 ⁇ p M, it is possible to set the p 1 and p M an arbitrary real number. In this way, it can be said that there are innumerable update patterns of the coefficient p.
  • the coefficient p in the simulated branching algorithm corresponds to the reciprocal of the temperature T in simulated annealing (SA). Therefore, the transition probability of the particles can be adjusted according to the coefficient p.
  • the coefficient p may be incremented at a constant rate according to the number of updates. Further, the coefficient p may be incremented by a different increment depending on the update timing. Further, the value of the coefficient p may be updated at some update timings, and the update of the value of the coefficient p may be skipped at the remaining update timings.
  • the update patterns of the coefficients p and a specify the region of the solution space to be searched and the granularity of the search. Therefore, it is necessary to adjust the update pattern of the coefficients p and a so that the global optimum solution can be easily found.
  • the simulated branch algorithm may be solved by the replica exchange method.
  • FIG. 8 conceptually shows the coefficient p used when performing the calculation by the replica exchange method.
  • the vertical axis of FIG. 8 shows the value of the coefficient p.
  • the horizontal axis of FIG. 8 is a virtual axis prepared for comparison with FIG. 7.
  • the replica exchange method when the replica exchange method is used, the number of combinations of a plurality of determined coefficients p (p (1) , p (2) , ..., P (M) ) is increased.
  • the 1st vector and the 2nd vector are updated.
  • the first vector and the second vector prepared for each coefficient p are referred to as replicas. Then, by exchanging the first vector and the second vector between the replicas, it is possible to realize the search for the solution in different regions in the solution space.
  • the replica exchange method it is possible to increase the probability of reaching the optimum solution of the problem or an approximate solution close to it by repeating the exchange of the first vector and the second vector.
  • a replica may be assigned to each processor of the calculation server. Further, a replica may be assigned to each calculation server in the information processing system. Also, the amount of computational resources allocated by the replica may be different. The types or combinations of computational resources allocated by the replica may be different. However, the method of allocating computational resources to each replica is not limited.
  • the calculations of the first vector, the second vector, and their functions in each replica can be performed in parallel (first level parallelization).
  • the update processing of the first vector and the second vector, each of which contains a plurality of elements can be performed in parallel (second level parallelization).
  • the particle size at which the first level of parallelization is performed and the particle size at which the second level of parallelization is performed in the computational resources may be different.
  • the particle size of the latter may be larger than the particle size of the former.
  • the first level of parallelization can be realized with the particle size of a calculation server and a virtual computer.
  • the second level of parallelization can be realized at the particle size of the processor.
  • ⁇ (p (t)) can be used as a (t).
  • the exchange probability between the replica m 1 and the replica m 1 + 1 is expressed by the following equation (9).
  • the data exchange circuit 51 can compare the value of the equation (9) with a pseudo-random number in the range of 0 or more and less than 1, and perform condition determination based on the magnitude relationship between the former value and the latter value. If the latter value is larger than the former value, the condition determination is positive, and the data exchange circuit 51 can exchange the first vector and the second vector between the replica m 1 and the replica m 1 + 1. it can. On the other hand, when the latter value is larger than the former value, the condition judgment is negative, and the exchange processing of the first vector and the second vector between the replica m 1 and the replica m 1 + 1 can be skipped. ..
  • each replica may be assigned to each arithmetic circuit.
  • the data exchange circuit executes the condition determination for each pair of arithmetic circuits, and when the conditional determination is positive, the exchange of the first vector and the second vector used for the calculation between the arithmetic circuits, or the first. It may be configured to perform at least one of the exchanges of one coefficient. Further, the data exchange circuit is configured to execute the condition determination based on the first vector (plurality of first variables), the second vector (plurality of second variables) and the first coefficient in the arithmetic circuit included in the pair. It may have been done. For example, the condition determination may be performed using the equation (9), or the condition determination may be performed using another equation.
  • Each replica may be assigned to each computer.
  • the second computer of the information processing system executes the condition determination for each pair of the first computers, and when the condition determination is affirmative, the first vector and the second computer used for the calculation between the first computers. It may be configured to perform at least one of the exchange of vectors or the exchange of first coefficients.
  • the first coefficient (for example, the coefficient p) used in each arithmetic circuit may have a geometric progression relationship with each other.
  • the process of solving is executed.
  • Each of the N variables X (m) i corresponds to the spin s i of the Ising model.
  • each of the N variables Y (m) i corresponds to momentum. It is assumed that the variables X (m) i and Y (m) i are all continuous variables.
  • the symplectic Euler method can be used to convert (10) above into a discrete recurrence formula, and the simulated branching algorithm can be calculated by the replica exchange method.
  • (11) below shows an example of a simulated branching algorithm after conversion to a recurrence formula.
  • the degree of parallelism of processing can be increased, and the optimum solution of the problem or an approximate solution close to it can be obtained in a relatively short time.
  • the method described above in which the combination of coefficients p is automatically determined, facilitates the implementation and execution of the simulated branching algorithm.
  • FIG. 9 shows an example of processing when the simulated branch algorithm is calculated by the replica exchange method. Hereinafter, the process will be described with reference to FIG.
  • the calculation server retrieves matrix J ij and vector h i corresponding to the problem, the counter t and the coefficient p (m), to initialize the a (m) (step S111).
  • the value of the coefficient p (m) can be set so as to form a geometric progression.
  • the counter t is set to, for example, 0.
  • the calculation server initializes the element X (m) i of the first vector and the element Y (m) i of the second vector, and substitutes 1 for m1 (step S112).
  • m1 is a parameter used to specify the replica m.
  • the initial values of the elements X (m) i of the first vector and the elements Y (m) i of the second vector can be set by pseudo-random numbers.
  • the calculation server updates the element X (m1) i of the first vector based on the corresponding element Y (m) i of the second vector (step S113). For example, in step S113, ⁇ t ⁇ D ⁇ Y (m1) i can be added to the element X (m1) i of the first vector. Then, the calculation server updates the element Y (m1) i of the second vector (step S114).
  • step S114 ⁇ t ⁇ [ (pDK ⁇ X (m1) i ⁇ X (m1) i ) ⁇ X (m1) i ] is added to the element Y (m1) i of the second vector. Can be done.
  • step S114 it is possible to add a - ⁇ t ⁇ c ⁇ h i ⁇ a - ⁇ t ⁇ c ⁇ ⁇ J ij ⁇ X (m1) i and into the variable Y (m1) i.
  • the calculation server can generate a pseudo-random number r of 0 or more and less than 1, and compare the value of r with the value of the above equation (9). When the value of the above equation (9) is larger than the value of the pseudo-random number r, it can be determined that the transition condition is satisfied. On the other hand, when the value of the pseudo-random number r is equal to or greater than the value of the above equation (9), it can be determined that the transition condition is not satisfied. Then, the calculation server updates the counter t (step S117). For example, in step S117, ⁇ t can be added to t.
  • step S118 it is determined whether or not the counter t is less than T (step S118). If the determination in step S118 is affirmative (YES in step S118), the calculation server re-executes the processes of steps S113 to S117. If the determination in step S118 is negative (NO in step S118), the compute server finds the spin s i , which is an element of the solution vector, based on the value of the element (variable) X (m) i of the first vector. (Step S119). For example, in step S119, the variable X (m) i is a positive value +1, may be converted respectively to -1 variable X (m) i is a negative value. In step S118, it may be determined whether or not the number of times the process of step S116 is executed exceeds the threshold value instead.
  • variable X (m) i in step S119 the representative value of the variable X (m1) i selected from any replica of 1 to M can be used.
  • the representative values are variable X (m1) selected from one of the replica of 1 ⁇ M may be a i, a variable X that is selected from at least one replica of 1 ⁇ M (m1) It may be the average value of i .
  • the information processing method may use a plurality of first computers and a second computer.
  • the plurality of first computers are configured to repeatedly update the first vector having the first variable as an element and the second vector having the second variable corresponding to the first variable as an element.
  • the second computer is configured to perform data exchange between the first computers.
  • the information processing method includes a step in which each first computer updates the first variable based on the corresponding second variable, and a second in which each first computer weights the first variable with a first coefficient.
  • the information processing method may use a plurality of first processors and a second processor.
  • the plurality of first processors are configured to repeatedly update the first vector having the first variable as an element and the second vector having the second variable corresponding to the first variable as an element.
  • the second processor is configured to perform the exchange of data between the first processors.
  • the information processing method includes a step in which each first processor updates the first variable based on the corresponding second variable, and a second in which each first processor weights the first variable with a first coefficient.
  • a step of adding to two variables a step of each first processor calculating a problem term using a plurality of first variables, a step of each first processor adding a problem term to a second variable, and a second.
  • the program according to the present embodiment may include the same steps as the information processing method. Further, the non-temporary computer-readable storage medium according to the present embodiment may store the program.
  • the information processing apparatus may vary the value of the coefficient p (m) based on the following equation (13).
  • Equation (13) is an example of a coefficient p (m) that can vary depending on the replica exchange timing or exchange determination timing.
  • the p (m) (t) of the equation (13) indicate that the value of the coefficient p (m) depends on the parameter t.
  • cos ( ⁇ ⁇ ⁇ ) of the formula (13), (- 1, + 1) to take a value in the range of, p (m) is .alpha.p 0 (m + 1) + .beta.p p as 0 (m) 0 (m + 1) and p 0 (m) each take a weighted value. Due to the periodicity of cos blunt, different values of p (m) are used for replicas with an even number m and replicas with an odd number m. In this way, each arithmetic circuit may periodically change the value of the first coefficient according to the number of updates of the first vector and the second vector.
  • the following equation (14) holds at the timing when the exchange determination between replicas is performed.
  • the relationship of the equation (14) is established at the exchange determination timing between the replicas, the first vector and the second vector are exchanged for the replica m 1 and the replica m 1 + 1 at each exchange determination timing. Since it is not necessary to perform the calculation shown in the equation (9) at each exchange determination timing, the consumption of computational resources can be suppressed. Further, since the determination process of whether or not to perform the exchange can be skipped at the exchange determination timing (step S116) in the flowchart of FIG. 9, the exchange determination timing of FIG. 9 is replaced with the exchange timing.
  • the user When using the coefficient p (m) shown in the equation (13), the user needs to determine the values of p 0 (m + 1) and p 0 (m) in advance. However, since it is not necessary to determine the update pattern of the coefficient p (m) and the definition of the coefficient a (m) as in the case of calculating the time evolution, it is possible to reduce the man-hours required to execute the simulated branch algorithm. Is. Further, by using the coefficient p (m) shown in the equation (13), the first vector and the second vector can be efficiently exchanged between the replicas, and the optimum solution of the problem or an approximation close thereto can be performed in a relatively short time. It is possible to obtain a solution.
  • FIG. 10 shows an example of processing when a calculation is performed using a variable coefficient p (m) .
  • p variable coefficient
  • the calculation server retrieves matrix J ij and vector h i corresponding to the problem, the counter t and the coefficient p (m), to initialize the a (m) (step S121).
  • the counter t is set to, for example, 0.
  • p (m) p 0 (m) is set.
  • step S125 If the determination in step S125 is affirmative (YES in step S125), the calculation server selects a random m1 from 0 to M every k times, and the values of X (m1) i and X (m1 + 1) i , The values of Y (m1) i and Y (m1 + 1) i are exchanged (step S126).
  • step S126 When the relationship of the equation (14) is established, it is not necessary to determine whether or not to exchange between the replicas, so step S126 corresponds to the exchange timing of the first vector and the second vector between the replicas. To do.
  • the calculation server updates the counter t (step S127). For example, in step S127, ⁇ t can be added to t.
  • step S128) it is determined whether or not the counter t is less than T (step S128). If the determination in step S128 is affirmative (YES in step S128), the calculation server re-executes the processes of steps S123 to S127. If the determination in step S128 is negative (NO in step S128), the compute server finds the spin s i , which is an element of the solution vector, based on the value of the element (variable) X (m) i of the first vector. (Step S129). For example, in step S129, the variable X (m) i is a positive value +1, may be converted respectively to -1 variable X (m) i is a negative value. It is the same as step S119 of FIG. 9 in that the method of selecting and calculating the representative value that can be used as the variable X (m) i in step S129 is not limited. Further, the process of step S129 may be executed by a server (computer) other than the calculation server.
  • the degree of parallelism of the processing is increased while suppressing the amount of calculation, and the problem occurs in a relatively short time. It is possible to obtain the optimum solution of or an approximate solution close to it.
  • the exchange processing of the first vector and the second vector between them (second exchange processing) can be executed alternately. This makes it possible to parallelize the exchange processing of the first vector and the second vector between the replicas.
  • FIG. 11 shows an example of processing in the case of performing replica exchange in parallel. Hereinafter, the process will be described with reference to FIG.
  • the calculation server retrieves matrix J ij and vector h i corresponding to the problem, the counter t and the coefficient p (m), to initialize the a (m) (step S131).
  • the method of initializing the coefficient p (m) in step S131 may be the same as that of step S121 described above.
  • the method of initializing the coefficient p (m) in step S131 may be the same as that of step S111 described above.
  • a combination of coefficients p (m) obtained on the empirical side may be used.
  • the counter t is set to, for example, 0.
  • the calculation server initializes the element X (m) i of the first vector and the element Y (m) i of the second vector, and substitutes 1 for m 1 (step S132).
  • m 1 is a parameter used to specify the replica m.
  • the initial values of the elements X (m) i of the first vector and the elements Y (m) i of the second vector can be set by pseudo-random numbers.
  • the calculation server updates the element X (m1) i of the first vector based on the corresponding element Y (m) i of the second vector (step S133).
  • step S133 ⁇ t ⁇ D ⁇ Y (m1) i can be added to the element X (m1) i of the first vector. Then, the calculation server updates the element Y (m1) i of the second vector (step S134). For example, in step S134, ⁇ t ⁇ [ (pdK ⁇ X (m1) i ⁇ X (m1) i ) ⁇ X (m1) i ] is added to the element Y (m1) i of the second vector. Can be done.
  • step S134 it is possible to further add the - ⁇ t ⁇ c ⁇ h i ⁇ a - ⁇ t ⁇ c ⁇ ⁇ J ij ⁇ X (m1) i in the element Y (m1) i of the second vector.
  • a m 1 is m 1 or odd an even number for each k times alternately selects determines a transition condition, if the transition condition is satisfied, X (m1) i and X
  • the value of (m1 + 1) i and the value of Y (m1) i and Y (m1 + 1) i are exchanged (step S135).
  • the calculation server can generate a pseudo-random number r of 0 or more and less than 1, and compare the value of r with the value of the above equation (9). When the value of the above equation (9) is larger than the value of the pseudo-random number r, it can be determined that the transition condition is satisfied. On the other hand, when the value of the pseudo-random number r is equal to or greater than the value of the above equation (9), it can be determined that the transition condition is not satisfied.
  • step S138 it is determined whether or not the counter t is less than T (step S138). If the determination in step S138 is affirmative (YES in step S138), the calculation server re-executes the processes of steps S133 to S137. If the determination in step S138 is negative (NO in step S138), the compute server finds the spin s i , which is an element of the solution vector, based on the value of the element (variable) X (m) i of the first vector. (Step S139). For example, in step S139, the variable X (m) i is a positive value +1, may be converted respectively to -1 variable X (m) i is a negative value. It is the same as in step S119 of FIG. 9 in that the method of selecting and calculating the representative value that can be used as the variable X (m) i in step S139 is not limited. Further, the process of step S139 may be executed by a server (information processing device) other than the calculation server.
  • the update processing of the elements of the first vector and the elements of the second vector is executed in parallel at least in a part of steps S132 to S135. May be (second level parallelization). Further, at least a part of the processing related to each replica may be parallelized (first level parallelization). For example, the processing of steps S132 to S135 relating to replicas different for each calculation server, virtual computer, or processor may be executed. In this case, the process of step S135 may be executed by a program running on any of the calculation servers, or may be executed by the control unit 13 of the management server 1 described above.
  • FIG. 12 conceptually shows the processing flow in the replica exchange method.
  • each square corresponds to a replica numbered 1-M.
  • Step S10 in FIG. 12 shows when the update processing of the first vector and the second vector is executed for each replica. Then, in step S11, it is determined whether or not the transition condition is satisfied for the even-numbered replica, and if the transition condition is satisfied, the exchange processing of the first vector and the second vector is executed. Further, in step S12, it is determined whether or not the transition condition is satisfied for the even-numbered replica, and if the transition condition is satisfied, the exchange processing of the first vector and the second vector is executed. As described above, the transition condition can be determined by generating a pseudo-random number r of 0 or more and less than 1, and comparing the value of r with the value of the above equation (9).
  • the data exchange circuit may be configured to select different pairs of arithmetic circuits depending on the timing and execute at least one of the exchange of the first vector and the second vector, or the exchange of the first coefficient. Further, the data exchange circuit alternately selects the pair of the first combination and the pair of the second combination in the plurality of arithmetic circuits, and exchanges the first vector and the second vector, or at least the exchange of the first coefficient. It may be configured to do either. For example, assume that there are replicas of numbers 1 to 6. It is assumed that the replicas of numbers 1 to 6 are assigned to any of the arithmetic circuits.
  • the first combination may include a pair of replicas of numbers 1 and 2, a pair of replicas of numbers 3 and 4, and a pair of replicas of numbers 5 and 6.
  • the second combination may include a pair of replicas of numbers 2 and 3, a pair of replicas of numbers 4 and 5, and a pair of replicas of numbers 6 and 1.
  • step S11 and the process of step S12 can be alternately executed depending on the timing.
  • the replica exchange process can be executed in parallel, and the optimum solution of the problem or an approximate solution close to the problem can be obtained in a shorter time.
  • FIG. 13 shows an example of an expression related to the replica exchange method.
  • Equation (E1) in FIG. 13 shows an example of the Hamiltonian used in each replica.
  • Equation (E2) shows an example of a canonical ensemble.
  • Equation (E3) shows a general detailed balance equation.
  • Equation (E4) shows detailed balance equations in canonical ensembles.
  • Equation (E5) is a summary of the terms of detailed balance equations in a canonical ensemble.
  • Equation (E6) shows the transition probability between replica m 1 and replica m 2 .
  • the graph of FIG. 14 shows the number of trials required to obtain the optimum solution in each of the case where the simulated branch algorithm is calculated by the replica exchange method and the case where the simulated branch algorithm is calculated by time evolution. It is a histogram showing an example of.
  • FIG. 14 shows the data when the Hamiltonian cycle problem of 48 nodes and 96 edges is solved by using the simulated branching algorithm.
  • the vertical axis of FIG. 14 shows the number of times that the optimum solution was obtained.
  • the calculation by the replica exchange method the calculation is performed using 10 replicas. Therefore, the number of trials in one calculation by the replica exchange method is counted as 10 times.
  • the replica exchange method it is determined whether or not the optimum solution is obtained every 50 cycles. Referring to the result of FIG. 14, when the replica exchange method is used, the optimum solution is obtained with a smaller number of trials than in the case of calculating the time evolution.
  • the information processing system can be made redundant by performing calculations using the replica exchange method. For example, even if the computer executing the calculation of the first vector and the second vector in one of the replicas stops abnormally, the computer executing the calculation of the first vector and the second vector in the other replica operates. If this is continued, the calculation process in the information processing system can be continued. In this case, the replica data used in the computer that stopped abnormally may be recovered based on the data in the other replicas. Even if the processor executing the calculation of the first vector and the second vector in any of the replicas stops abnormally, the calculation can be continued in the other processor.
  • Equation (15) corresponds to the energy of the Ising model including many-body interactions.
  • both QUABO and HOBO are a kind of unconstrained polynomial binary variable optimization (PUBO: Polynomial Unconstrained Binary Optimization). That is, among PUBOs, the combinatorial optimization problem having a quadratic objective function is QUABO. Further, among PUBOs, the combinatorial optimization problem having a third-order or higher objective function can be said to be HOBO.
  • PUBO Polynomial Unconstrained Binary Optimization
  • the Hamiltonian H in the above equation (5) may be replaced with the Hamiltonian H in the following equation (16).
  • the problem term shown in the following equation (17) is derived from the equation (16).
  • the problem item z i of (17) takes the form of partially differentiating the second equation of (16) with respect to any variable x i (element of the first vector).
  • the variable x i to be partially differentiated depends on the index i.
  • the index i of the variable x i corresponds to the index that specifies the element of the first vector and the element of the second vector.
  • the problem items shown above are only examples of problem items that can be used by the information processing apparatus according to the present embodiment. Therefore, the form of the problem argument used in the calculation may be different from these.
  • additional processing may be performed when updating the elements of the first vector in order to reduce calculation errors.
  • the absolute value of the element X (n) i of the first vector is larger than 1 by the update replaces the value of the element X (n) i of the first vector to sgn (X (n) i) . That is, when X (n) i > 1 due to the update, the value of the variable X (n) i is set to 1. Further, when X (n) i ⁇ -1 due to the update, the value of the variable X (n) i is set to -1. This makes it possible to approximate the spin s i with higher accuracy using the variable X (n) i .
  • X (n) i> 1 by the update may be multiplied by a coefficient rf variable Y (n) i corresponding to the variable X (n) i.
  • a coefficient rf of -1 ⁇ r ⁇ 0 is used, the above wall becomes a wall having a reflection coefficient rf.
  • the arithmetic circuit sets the second variable corresponding to the first variable whose value is smaller than the first value, or the second variable corresponding to the first variable whose value is larger than the second value, to the original first variable. It may be configured to update to a value obtained by multiplying two variables by a second coefficient. For example, the arithmetic circuit sets the second variable corresponding to the first variable having a value smaller than -1 or the second variable corresponding to the first variable having a value greater than 1 as the original second variable with a second coefficient. It may be configured to update to a value multiplied by.
  • the second coefficient corresponds to the above-mentioned coefficient rf.
  • the arithmetic circuit when it becomes X (n) i> 1 by updating the value of the variable Y (n) i corresponding to the variable X (n) i may be set to a pseudo-random number.
  • a pseudo-random number For example, random numbers in the range [-0.1, 0.1] can be used. That is, the arithmetic circuit converts the value of the second variable corresponding to the first variable whose value is smaller than the first value, or the value of the second variable corresponding to the first variable whose value is larger than the second value into a pseudo random number. It may be configured to be set.
  • a continuous variable X is used in the problem term instead of a discrete variable. Therefore, there is a possibility that an error may occur with the discrete variable used in the original combinatorial optimization problem.
  • the value sgn (X) obtained by converting the continuous variable X with a sign function can be used instead of the continuous variable X in the calculation of the problem term as shown in (20) below.
  • sgn (X) corresponds to spin s.
  • the product of the spins appearing in the problem term always takes a value of either -1 or 1, so when dealing with the HOMO problem having a higher-order objective function, an error occurs due to the product operation. Can be prevented.
  • a spin vector can be obtained by converting each element of the first vector with a sign function.
  • a PC cluster is a system that connects a plurality of computers and realizes computing performance that cannot be obtained by one computer.
  • the information processing system 100 shown in FIG. 1 includes a plurality of calculation servers and processors, and can be used as a PC cluster.
  • MPI Message Passing Interface
  • MPI can be used to implement the control program 14E of the management server 1, the calculation program 34B and the control program 34C of each calculation server.
  • the calculation of L variables among the variables x i included in the first vector (x 1 , x 2 , ..., X N ) for each processor. Can be done.
  • m (j-1) L + 1, (j-1) L + 2, ..., jL ⁇ and ⁇ y.
  • m ( j-1) L + 1, (j-1) L + 2, ⁇ , jL ⁇ calculated tensor shown in (21) as follows required to J of (n ) Shall be stored in a storage area accessible to processor # j (eg, register, cache, memory, etc.).
  • each processor calculates a constant variable of the first vector and the second vector.
  • the number of variables in the first vector and the second vector to be calculated may differ depending on the processor. For example, when there is a performance difference depending on the processor implemented in the calculation server, the number of variables to be calculated can be determined according to the performance of the processor.
  • the values of all the components of the first vector are required.
  • the conversion to a binary variable can be performed, for example, by using the sign function sgn (). Therefore, using the Allgather function, the values of all the components of the first vector (x 1 , x 2 , ..., X N ) can be shared by the Q processors.
  • the second vector (y 1 , y 2 , ..., Y N )
  • tensor J (n) it is not essential to share values between processors.
  • Data sharing between processors can be realized, for example, by using interprocess communication or storing data in shared memory.
  • Processor # j calculates the value of the problem item ⁇ z m
  • m (j-1) L + 1, (j-1) L + 2, ..., JL ⁇ .
  • the processor #j is calculated problem claim ⁇ z m
  • m ( j-1) L + 1, based on the value of (j-1) L + 2 , ⁇ , jL ⁇ , variables ⁇ y m
  • m Update (j-1) L + 1, (j-1) L + 2, ..., JL ⁇ .
  • the tensor J (n) and the vector (x 1 , x 2 , ... , X N ) requires a product-sum operation, including the calculation of the product.
  • FIG. 15 schematically shows an example of a multiprocessor configuration.
  • the plurality of calculation nodes in FIG. 15 correspond to, for example, a plurality of calculation servers of the information processing system 100.
  • the high-speed link in FIG. 15 corresponds to, for example, an interconnect between cables 4a to 4c of the information processing system 100 and a calculation server formed by the switch 5.
  • the shared memory of FIG. 15 corresponds to, for example, the shared memory 32.
  • the processor of FIG. 15 corresponds to, for example, the processors 33A to 33D of each calculation server.
  • FIG. 15 shows the data arranged in each component and the data transferred between the components.
  • the values of the variables x i and y i are calculated. Further, between the processor and the shared memory, the variable x i is transferred.
  • L of the first vector (x 1 , x 2 , ..., X N ) and the second vector (y 1 , y 2 , ..., Y N ) are stored.
  • the variable and a part of the tensor J (n) are saved.
  • the first vector (x 1 , x 2 , ..., X N ) is transferred.
  • all the elements of the first vector (x 1 , x 2 , ..., X N ) are required to update the variable y i in each processor.
  • one or more calculation nodes can be assigned for the calculation of the first vector and the second vector in each replica.
  • the data of the first vector and the second vector are transferred between the calculation nodes via the high-speed link.
  • the control unit 13 of the management server 1 can control the calculation node.
  • the arrangement and transfer of data shown in FIG. 15 is only an example. The method of arranging data, the method of transferring data, and the method of realizing parallelization in a PC cluster are not particularly limited.
  • the simulated branch algorithm may be calculated using the GPU (Graphics Processing Unit).
  • FIG. 16 schematically shows an example of a configuration using a GPU.
  • FIG. 16 shows a plurality of GPUs connected to each other by a high-speed link.
  • Each GPU has a plurality of cores that can access the shared memory.
  • a plurality of GPUs are connected via a high-speed link to form a GPU cluster.
  • the high speed link corresponds to an interconnect between the compute servers formed by the cables 4a-4c and the switch 5.
  • parallel calculations can be executed even when one GPU is used. That is, each GPU in FIG. 16 can execute the calculation corresponding to each calculation node in FIG. That is, the processor of the information processing device (calculation server) may be the core of the Graphics Processing Unit (GPU).
  • the processor of the information processing device may be the core of the Graphics Processing Unit (GPU).
  • the variables x i and y i , and the tensor J (n) are defined as device variables.
  • the GPU can calculate the product of the tensor J (n) required to update the variable y i and the first vector (x 1 , x 2 , ..., X N ) in parallel by the matrix vector product function. it can.
  • the product of the tensor and the vector can be obtained by repeatedly executing the product operation of the matrix and the vector.
  • the calculation of the first vector (x 1 , x 2 , ..., X N ) and the part of the second vector (y 1 , y 2 , ..., y N ) other than the multiply-accumulate operation. is, i-th element in each thread (x i, y i) to execute the update processing, it is possible to realize parallel processing.
  • FIG. 17 shows an example of the overall processing executed to solve the combinatorial optimization problem. Hereinafter, the process will be described with reference to FIG.
  • the combinatorial optimization problem is formulated (step S201). Then, the formulated combinatorial optimization problem is converted into an Ising problem (Ising model format) (step S202). Next, the solution of the Ising problem is calculated by the Ising machine (information processing device) (step S203). Then, the calculated solution is verified (step S204). For example, in step S204, it is confirmed whether or not the constraint condition is satisfied. Further, in step S204, the value of the energy function (Hamiltonian) may be referred to, and it may be confirmed whether or not the obtained solution is the optimum solution or an approximate solution close to it.
  • the energy function Haltonian
  • step S205 it is determined whether or not to recalculate according to at least one of the verification result in step S204 and the number of calculations. If it is determined that the recalculation is to be performed (YES in step S205), the processes of steps S203 and S204 are executed again. On the other hand, when it is determined not to recalculate (NO in step S205), a solution is selected (step S206). For example, in step S206, selection can be made based on at least either the satisfaction of the constraints or the value of the energy function. If a plurality of solutions have not been calculated, the process of step S206 may be skipped. Finally, the selected solution is converted into the solution of the combinatorial optimization problem, and the solution of the combinatorial optimization problem is output (step S207).
  • the present invention is not limited to the above embodiment as it is, and at the implementation stage, the components can be modified and embodied within a range that does not deviate from the gist thereof.
  • various inventions can be formed by an appropriate combination of the plurality of components disclosed in the above-described embodiment. For example, some components may be removed from all the components shown in the embodiments. In addition, components across different embodiments may be combined as appropriate.
  • Management server 2 Network 3a, 3b, 3c Calculation server 4a, 4b, 4c Cable 5 Switch 6 Client terminal 10 Processor 11 Management unit 12 Conversion unit 13 Control unit 14 Storage unit 14A Problem data 14B Calculation data 14C Management program 14D Conversion program 14E , 34C Control program 15, 31 Communication circuit 16 Input circuit 17 Output circuit 18 Operation device 19 Display device 20 Bus 32 Shared memory 33A, 33B, 33C, 33D Processor 34 Storage 34A Calculation data 34B Calculation program 35 Host bus adapter

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Operations Research (AREA)
  • Computational Linguistics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

［課題］組合せ最適化問題の解を実用的な時間内で計算する情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラムを提供する。 ［解決手段］　本発明の実施形態としての情報処理装置は、それぞれが第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている、複数の演算回路と、データ交換回路とを備えている。それぞれの前記演算回路は、前記第１変数を対応する第２変数に基づいて更新し、前記第１変数を第１係数で重み付けし対応する前記第２変数に加算し、複数の前記第１変数を用いて問題項を計算し、前記問題項を前記第２変数に加算するように構成されており、それぞれの前記演算回路において前記第１係数に異なる値が設定されており、前記データ交換回路は、前記演算回路間で前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されている。

Description

情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム
　本発明の実施形態は、情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラムに関する。
　組合せ最適化問題とは、複数の組合せの中から目的に最も適した組合せを選ぶ問題である。組合せ最適化問題は、数学的には、「目的関数」と呼ばれる、複数の離散変数を有する関数を最大化させる問題、または、当該関数を最小化させる問題に帰着される。組合せ最適化問題は、金融、物流、交通、設計、製造、生命科学など各種の分野において普遍的な問題であるが、組合せ数が問題サイズの指数関数のオーダーで増える、いわゆる「組合せ爆発」のため、必ず最適解を求めることができるとは限らない。また、最適解に近い近似解を得ることすら難しい場合が多い。
　各分野における問題を解決し、社会のイノベーションおよび科学技術の進歩を促進するために、組合せ最適化問題の解を実用的な時間内で計算する技術の開発が求められている。
特開２０１７－７３１０６号公報
H. Goto, K. Tatsumura, A. R. Dixon, Sci. Adv. 5, eaav2372 (2019). H. Goto, Sci. Rep. 6, 21686 (2016). 土屋、西山、辻田：分岐特性を用いた組合せ最適化問題の近似解法URL:http://www.ynl.t.u-tokyo.ac.jp/project/RobotBrainCREST/publications/pdf/tsuchiya/4_01.pdf 土屋、西山、辻田：決定論的アニーリングアルゴリズムの解析URL:http://www.ynl.t.u-tokyo.ac.jp/project/RobotBrainCREST/publications/pdf/tsuchiya/4_02.pdf
　本発明の実施形態は、組合せ最適化問題の解を実用的な時間内で計算する情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラムを提供する。
　本発明の実施形態としての情報処理装置は、それぞれが第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている、複数の演算回路と、データ交換回路とを備えている。それぞれの前記演算回路は、前記第１変数を対応する第２変数に基づいて更新し、前記第１変数を第１係数で重み付けし対応する前記第２変数に加算し、複数の前記第１変数を用いて問題項を計算し、前記問題項を前記第２変数に加算するように構成されており、それぞれの前記演算回路において前記第１係数に異なる値が設定されており、前記データ交換回路は、前記演算回路間で前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されている。
情報処理システムの構成例を示した図。 管理サーバの構成例を示したブロック図。 管理サーバの記憶部に保存されるデータの例を示す図。 計算サーバの構成例を示したブロック図。 計算サーバのストレージに保存されるデータの例を示す図。 時間発展によってシミュレーテッド分岐アルゴリズムの解を計算する場合における処理の例を示したフローチャート。 時間発展の計算を行う場合における係数ｐの更新パターンの例を示したグラフ。 レプリカ交換法によって計算を行う場合に使われる係数ｐを概念的に示した図。 レプリカ交換法によってシミュレーテッド分岐アルゴリズムの計算を行ったときにおける処理の例を示したフローチャート。 変動する係数ｐ（ｍ）を使って計算を行ったときにおける処理の例を示したフローチャート。 並列的にレプリカ交換を行う場合における処理の例を示したフローチャート。 レプリカ交換法における処理フローを概念的に示した図。 レプリカ交換法に関連する式の例を示した図。 最適解を得られるまでに必要な試行回数の例を示したヒストグラム。 マルチプロセッサ構成の例を概略的に示した図。 ＧＰＵを使った構成の例を概略的に示した図。 組合せ最適化問題を解くために実行される全体的な処理の例を示したフローチャート。
　以下、図面を参照しながら、本発明の実施形態について説明する。また、図面において同一の構成要素は、同じ番号を付し、説明は、適宜省略する。
　図１は、情報処理システム１００の構成例を示したブロック図である。図１の情報処理システム１００は、管理サーバ１と、ネットワーク２と、計算サーバ（情報処理装置）３ａ～３ｃと、ケーブル４ａ～４ｃと、スイッチ５と、記憶装置７を備えている。また、図１には、情報処理システム１００と通信可能なクライアント端末６が示されている。管理サーバ１、計算サーバ３ａ～３ｃ、クライアント端末６および記憶装置７は、ネットワーク２を介して互いにデータ通信をすることができる。例えば、計算サーバ３ａ～３ｃは、記憶装置７にデータを保存したり、記憶装置７よりデータを読み出したりすることができる。ネットワーク２は、例えば、複数のコンピュータネットワークが相互に接続されたインターネットである。ネットワーク２は、通信媒体として有線、無線、または、これらの組合せを用いることができる。また、ネットワーク２で使われる通信プロトコルの例としては、ＴＣＰ／ＩＰがあるが、通信プロトコルの種類については特に問わない。
　また、計算サーバ３ａ～３ｃは、それぞれケーブル４ａ～４ｃを介してスイッチ５に接続されている。ケーブル４ａ～４ｃおよびスイッチ５は、計算サーバ間のインターコネクトを形成している。計算サーバ３ａ～３ｃは、当該インターコネクトを介して互いにデータ通信をすることも可能である。スイッチ５は、例えば、Ｉｎｆｉｎｉｂａｎｄのスイッチである。ケーブル４ａ～４ｃは、例えば、Ｉｎｆｉｎｉｂａｎｄのケーブルである。ただし、Ｉｎｆｉｎｉｂａｎｄのスイッチ／ケーブルの代わりに、有線ＬＡＮのスイッチ／ケーブルを使ってもよい。ケーブル４ａ～４ｃおよびスイッチ５で使われる通信規格および通信プロトコルについては、特に問わない。クライアント端末６の例としては、ノートＰＣ、デスクトップＰＣ、スマートフォン、タブレット、車載端末などが挙げられる。
　組合せ最適化問題の求解では、処理の並列化および／または処理の分散化を行うことができる。したがって、計算サーバ３ａ～３ｃおよび／または計算サーバ３ａ～３ｃのプロセッサは、計算処理の一部のステップを分担して実行してもよいし、異なる変数について同様の計算処理を並列的に実行してもよい。管理サーバ１は、例えば、ユーザによって入力された組合せ最適化問題を各計算サーバが処理可能な形式に変換し、計算サーバを制御する。そして、管理サーバ１は、各計算サーバから計算結果を取得し、集約した計算結果を組合せ最適化問題の解に変換する。こうして、ユーザは、組合せ最適化問題の解を得ることができる。組合せ最適化問題の解は、最適解と、最適解に近い近似解とを含むものとする。
　図１には、３台の計算サーバが示されている。ただし、情報処理システムに含まれる計算サーバの台数を限定するものではない。また、組合せ最適化問題の求解に使われる計算サーバの台数についても特に問わない。例えば、情報処理システムに含まれる計算サーバは１台であってもよい。また、情報処理システムに含まれる複数の計算サーバのうち、いずれかの計算サーバを使って組合せ最適化問題の求解を行ってもよい。また、情報処理システムに、数百台以上の計算サーバが含まれていてもよい。計算サーバは、データセンターに設置されたサーバであってもよいし、オフィスに設置されたデスクトップＰＣであってもよい。また、計算サーバは異なるロケーションに設置された複数の種類のコンピュータであってもよい。計算サーバとして使われる情報処理装置の種類については特に問わない。例えば、計算サーバは、汎用的なコンピュータであってもよいし、専用の電子回路または、これらの組合せであってもよい。
　図２は、管理サーバ１の構成例を示したブロック図である。図２の管理サーバ１は、例えば、中央演算処理装置（ＣＰＵ）とメモリとを含むコンピュータである。管理サーバ１は、プロセッサ１０と、記憶部１４と、通信回路１５と、入力回路１６と、出力回路１７とを備えている。プロセッサ１０、記憶部１４、通信回路１５、入力回路１６および出力回路１７は、互いにバス２０を介して接続されているものとする。プロセッサ１０は、内部の構成要素として、管理部１１と、変換部１２と、制御部１３とを含んでいる。
　プロセッサ１０は、演算を実行し、管理サーバ１の制御を行う電子回路である。プロセッサ１０として、例えば、ＣＰＵ、マイクロプロセッサ、ＡＳＩＣ、ＦＰＧＡ、ＰＬＤまたはこれらの組合せを用いることができる。管理部１１は、ユーザのクライアント端末６を介して管理サーバ１の操作を行うためのインタフェースを提供する。管理部１１が提供するインタフェースの例としては、ＡＰＩ、ＣＬＩ、ウェブページなどが挙げられる。例えば、ユーザは、管理部１１を介して組合せ最適化問題の情報の入力を行ったり、計算された組合せ最適化問題の解の閲覧および／またはダウンロードを行ったりすることができる。変換部１２は、組合せ最適化問題を各計算サーバが処理可能な形式に変換する。制御部１３は、各計算サーバに制御指令を送信する。制御部１３が各計算サーバから計算結果を取得した後、変換部１２は、複数の計算結果を集約し、組合せ最適化問題の解に変換する。また、制御部１３は、各計算サーバまたは各サーバ内のプロセッサが実行する処理内容を指定してもよい。
　記憶部１４は、管理サーバ１のプログラム、プログラムの実行に必要なデータおよびプログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、ＯＳとアプリケーションの両方を含むものとする。記憶部１４は、揮発性メモリ、不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリの例としては、ＤＲＡＭ、ＳＲＡＭなどがある。不揮発性メモリの例としては、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭ、またはＭＲＡＭが挙げられる。また、記憶部１４として、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置を使ってもよい。
　通信回路１５は、ネットワーク２に接続された各装置との間でデータの送受信を行う。通信回路１５は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ　Ｉｎｔｅｒｆａｃｅ　Ｃａｒｄ）である。ただし、通信回路１５は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。入力回路１６は、管理サーバ１へのデータ入力を実現する。入力回路１６は、外部ポートとして、例えば、ＵＳＢ、ＰＣＩ－Ｅｘｐｒｅｓｓなどを備えているものとする。図２の例では、操作装置１８が入力回路１６に接続されている。操作装置１８は、管理サーバ１に情報を入力するための装置である。操作装置１８は、例えば、キーボード、マウス、タッチパネル、音声認識装置などであるが、これに限られない。出力回路１７は、管理サーバ１からのデータ出力を実現する。出力回路１７は、外部ポートとしてＨＤＭＩ、ＤｉｓｐｌａｙＰｏｒｔなどを備えているものとする。図２の例では、表示装置１９が出力回路１７に接続されている。表示装置１９の例としては、ＬＣＤ（液晶ディスプレイ）、有機ＥＬ（有機エレクトロルミネッセンス）ディスプレイ、またはプロジェクタがあるが、これに限られない。
　管理サーバ１の管理者は、操作装置１８および表示装置１９を使って、管理サーバ１のメンテナンスを行うことができる。なお、操作装置１８および表示装置１９は、管理サーバ１に組み込まれたものであってもよい。また、管理サーバ１に必ず操作装置１８および表示装置１９が接続されていなくてもよい。例えば、管理者は、ネットワーク２と通信可能なクライアント端末を用いて管理サーバ１のメンテナンスを行ってもよい。
　図３は、管理サーバ１の記憶部１４に保存されるデータの例を示している。図３の記憶部１４には、問題データ１４Ａと、計算データ１４Ｂと、管理プログラム１４Ｃと、変換プログラム１４Ｄと、制御プログラム１４Ｅとが保存されている。例えば、問題データ１４Ａは、組合せ最適化問題のデータを含む。例えば、計算データ１４Ｂは、各計算サーバから収集された計算結果を含む。例えば、管理プログラム１４Ｃは、上述の管理部１１の機能を実現するプログラムである。例えば、変換プログラム１４Ｄは、上述の変換部１２の機能を実現するプログラムである。例えば、制御プログラム１４Ｅは、上述の制御部１３の機能を実現するプログラムである。
　図４は、計算サーバの構成例を示したブロックである。図４には、例示的に計算サーバ３ａの構成が示されている。他の計算サーバは、計算サーバ３ａと同様の構成であってもよいし、計算サーバ３ａと異なる構成であってもよい。計算サーバ３ａは、例えば、第１ベクトルおよび第２ベクトルの計算を単独で、または、他の計算サーバと分担して実行する情報処理装置である。また、少なくともいずれかの計算サーバは、第１ベクトルの要素間の問題項を計算してもよい。ここで、問題項は、後述するイジングモデルに基づいて計算されるものであってもよい。また、問題項は、後述する多体相互作用を含むものであってもよい。
　例えば、第１ベクトルは、変数ｘｉ（ｉ＝１、２、・・・、Ｎ）または変数Ｘ（ｍ） ｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルである。例えば、第２ベクトルは、変数ｙｉ（ｉ＝１、２、・・・、Ｎ）または変数Ｙ（ｍ）ｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルである。ここで、ｍは、１～Ｍの整数のレプリカ番号を示している。レプリカの詳細については、後述する。
　計算サーバ３ａは、例えば、通信回路３１と、共有メモリ３２と、プロセッサ３３Ａ～３３Ｄと、ストレージ３４と、ホストバスアダプタ３５とを備えている。通信回路３１、共有メモリ３２、プロセッサ３３Ａ～３３Ｄ、ストレージ３４、ホストバスアダプタ３５は、バス３６を介して互いに接続されているものとする。
　通信回路３１は、ネットワーク２に接続された各装置との間でデータの送受信を行う。通信回路３１は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ　Ｉｎｔｅｒｆａｃｅ　Ｃａｒｄ）である。ただし、通信回路３１は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。共有メモリ３２は、プロセッサ３３Ａ～３３Ｄからアクセス可能なメモリである。共有メモリ３２の例としては、ＤＲＡＭ、ＳＲＡＭなどの揮発性メモリが挙げられる。ただし、共有メモリ３２として、不揮発性メモリなどその他の種類のメモリが使われてもよい。プロセッサ３３Ａ～３３Ｄは、共有メモリ３２を介してデータの共有を行うことができる。なお、必ず計算サーバ３ａのすべてのメモリが共有メモリとして構成されていなくてもよい。例えば、計算サーバ３ａの一部のメモリは、いずれかのプロセッサのみからアクセスできるローカルメモリとして構成されていてもよい。
　プロセッサ３３Ａ～３３Ｄは、計算処理を実行する電子回路である。プロセッサは、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）、ＧＰＵ（Ｇｒａｐｈｉｃｓ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、またはＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）のいずれであってもよいし、これらの組合せであってもよい。また、プロセッサは、ＣＰＵコアまたはＣＰＵスレッドであってもよい。プロセッサがＣＰＵである場合、計算サーバ３ａが備えるソケット数については、特に問わない。また、プロセッサは、ＰＣＩ　ｅｘｐｒｅｓｓなどのバスを介して計算サーバ３ａのその他の構成要素に接続されていてもよい。
　図４の例では、計算サーバが４つのプロセッサを備えている。ただし、１台の計算サーバが備えているプロセッサの数はこれとは異なっていてもよい。例えば、計算サーバに実装されているプロセッサの数および／または種類が異なっていてもよい。
　情報処理装置は、例えば、第１変数ｘｉ（ｉ＝１、２、・・・、Ｎ）を要素とする第１ベクトルおよび第１変数に対応する第２変数ｙｉ（ｉ＝１、２、・・・、Ｎ）を要素とする第２ベクトルを繰り返し更新するように構成されている。情報処理装置の記憶部は、第１ベクトルの要素である第１変数および第２ベクトルの要素である第２変数を記憶するように構成されていてもよい。
　例えば、複数の演算回路５０は、それぞれが第１変数を要素とする第１ベクトルおよび第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている。例えば、それぞれの演算回路５０は、第１変数を対応する第２変数に基づいて更新し、第１変数を第１係数で重み付けし対応する第２変数に加算し、複数の第１変数を用いて問題項を計算し、問題項を第２変数に加算するように構成されている。また、それぞれの演算回路５０において第１係数に異なる値が設定されていてもよい。例えば、データ交換回路５１は、演算回路５０間で第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するように構成されている。演算回路が計算する問題項は、イジングモデルに基づいていてもよい。また、演算回路が計算する問題項は、多体相互作用を含んでいてもよい。問題項の詳細については、後述する。
　図４の例では、プロセッサ３３Ａ～３３Ｃが複数の演算回路に相当しており、プロセッサ３３Ｄがデータ交換回路に相当している。ただし、図４に示した演算回路／データ交換回路とプロセッサとの対応関係は、一例にしかすぎない。したがって、演算回路／データ交換回路とプロセッサとの対応関係はこれとは異なっていてもよい。また、演算回路／データ交換回路に割り当てられるプロセッサ数については、特に限定しない。後述するように、同一のプロセッサが演算回路およびデータ交換回路の役割を兼ね備えていてもよい。計算サーバに複数の種類のプロセッサ（例えば、ＣＰＵ、ＧＰＵ、ＦＰＧＡ）が実装されている場合には、異なる種類のプロセッサを演算回路およびデータ交換回路に割り当ててもよい。
　ここでは、プロセッサ単位に処理内容の割り当てが行われる場合を例に説明した。ただし、処理内容の割り当てが行われる計算資源の単位を限定するものではない。例えば、計算機単位で処理内容の割り当てを行ってもよい。この場合、情報処理システムは、複数の第１計算機と、第２計算機とを備えていてもよい。例えば、第１計算機は、それぞれが第１変数を要素とする第１ベクトルおよび第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている。例えば、第２計算機は、第１計算機間におけるデータの交換を実行するように構成されている。第１計算機と、第２計算機は、例えば、上述の計算サーバに相当する。
　それぞれの第１計算機は、第１変数を対応する第２変数に基づいて更新し、第１変数を第１係数で重み付けし対応する第２変数に加算し、複数の第１変数を用いて問題項を計算し、問題項を第２変数に加算するように構成されていてもよい。第２計算機は、第１計算機間で第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するように構成されていてもよい。また、それぞれの第１計算機において第１係数に異なる値が設定されていてもよい。
　以下では、再び図４を参照し、計算サーバの構成要素を説明する。
　ストレージ３４は、計算サーバ３ａのプログラム、プログラムの実行に必要なデータ、およびプログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、ＯＳとアプリケーションの両方を含むものとする。ストレージ３４は、揮発性メモリ、不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリの例としては、ＤＲＡＭ、またはＳＲＡＭなどがある。不揮発性メモリの例としては、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭ、またはＭＲＡＭが挙げられる。また、ストレージ３４として、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置が使われてもよい。
　ホストバスアダプタ３５は、計算サーバ間のデータ通信を実現する。ホストバスアダプタ３５は、ケーブル４ａを介してスイッチ５に接続されている。ホストバスアダプタ３５は、例えば、ＨＣＡ（Ｈｏｓｔ　Ｃｈａｎｎｅｌ　Ａｄａｐｔｏｒ）である。ホストバスアダプタ３５、ケーブル４ａ、およびスイッチ５を使って高スループットを実現可能なインターコネクトを形成することにより、並列的な計算処理の速度を向上させることができる。
　図５は、計算サーバのストレージに保存されるデータの例を示している。図５のストレージ３４には、計算データ３４Ａと、計算プログラム３４Ｂと、制御プログラム３４Ｃとが保存されている。計算データ３４Ａは、計算サーバ３ａの計算途中のデータまたは計算結果を含んでいる。なお、計算データ３４Ａの少なくとも一部は、共有メモリ３２、プロセッサのキャッシュ、またはプロセッサのレジスタなど、異なる記憶階層に保存されていてもよい。計算プログラム３４Ｂは、所定のアルゴリズムに基づき、各プロセッサにおける計算処理および、共有メモリ３２およびストレージ３４へのデータの保存処理を実現するプログラムである。制御プログラム３４Ｃは、管理サーバ１の制御部１３から送信された指令に基づき、計算サーバ３ａを制御し、計算サーバ３ａの計算結果を管理サーバ１に送信するプログラムである。
　次に、組合せ最適化問題の求解に関連する技術について説明する。組合せ最適化問題を解くために使われる情報処理装置の一例として、イジングマシンが挙げられる。イジングマシンとは、イジングモデルの基底状態のエネルギーを計算する情報処理装置のことをいう。これまで、イジングモデルは、主に強磁性体や相転移現象のモデルとして使われることが多かった。しかし、近年、イジングモデルは、組合せ最適化問題を解くためのモデルとしての利用が増えている。下記の式（１）は、イジングモデルのエネルギーを示している。
Figure JPOXMLDOC01-appb-M000001
 ここで、ｓｉ、ｓｊはスピンである、スピンは、＋１または－１のいずれかの値をとる２値変数である。Ｎは、スピンの数である。ｈｉは、各スピンに作用する局所磁場である。Ｊは、スピン間における結合係数の行列である。行列Ｊは、対角成分が０である実対称行列となっている。したがって、Ｊｉｊは行列Ｊのｉ行ｊ列の要素を示している。なお、式（１）のイジングモデルは、スピンについての２次式となっているが、スピンの３次以上の項を含む拡張されたイジングモデル（多体相互作用を有するイジングモデル）を使ってもよい。
　式（１）のイジングモデルを使うと、エネルギーＥＩｓｉｎｇを目的関数とし、エネルギーＥＩｓｉｎｇを可能な限り小さくする解を計算することができる。イジングモデルの解は、スピンのベクトル（ｓ１、ｓ２、・・・、ｓＮ）の形式で表される。特に、エネルギーＥＩｓｉｎｇが最小値となるベクトル（ｓ１、ｓ２、・・・、ｓＮ）は、最適解とよばれる。ただし、計算されるイジングモデルの解は、必ず厳密な最適解でなくてもよい。以降では、イジングモデルを使ってエネルギーＥＩｓｉｎｇが可能な限り小さくなる近似解（すなわち、目的関数の値が可能な限り最適値に近くなる近似解）を求める問題をイジング問題とよぶものとする。
　式（１）のスピンｓｉは２値変数であるため、式（１＋ｓｉ）／２を使うことにより、組合せ最適化問題で使われる離散変数（ビット）との変換を容易に行うことができる。したがって、組合せ最適化問題をイジング問題に変換し、イジングマシンに計算を行わせることにより、組合せ最適化問題の解を求めることが可能である。０または１のいずれかの値をとる離散変数（ビット）を変数とする２次の目的関数を最小化する解を求める問題は、ＱＵＢＯ（Ｑｕａｄｒａｔｉｃ　Ｕｎｃｏｎｓｔｒａｉｎｅｄ　Ｂｉｎａｒｙ　Ｏｐｔｉｍｉｚａｔｉｏｎ、制約なし２値変数２次最適化）問題とよばれる。式（１）で表されるイジング問題は、ＱＵＢＯ問題と等価であるといえる。
　例えば、量子アニーラ、コヒーレントイジングマシン、または量子分岐マシンなどがイジングマシンのハードウェア実装として提案されている。量子アニーラは、超伝導回路を使って量子アニーリングを実現する。コヒーレントイジングマシンは、光パラメトリック発振器で形成されたネットワークの発振現象を利用する。量子分岐マシンは、カー効果を有するパラメトリック発振器のネットワークにおける量子力学的な分岐現象を利用する。これらのハードウェア実装は、計算時間の大幅な短縮を実現する可能性がある一方、大規模化や安定的な運用が難しいという課題もある。
　そこで、広く普及しているデジタルコンピュータを使ってイジング問題の求解を行うことも可能である。デジタルコンピュータは、上述の物理的現象を使ったハードウェア実装と比べ、大規模化と安定運用が容易である。デジタルコンピュータでイジング問題の求解を行うためのアルゴリズムの一例として、シミュレーテッドアニーリング（ＳＡ）が挙げられる。シミュレーテッドアニーリングをより高速に実行する技術の開発が行われている。ただし、一般のシミュレーテッドアニーリングはそれぞれの変数が逐次更新される逐次更新アルゴリズムであるため、並列化による計算処理の高速化は難しい。
　上述の課題を踏まえ、デジタルコンピュータにおける並列的な計算によって、規模の大きい組合せ最適化問題の求解を高速に行うことが可能なシミュレーテッド分岐アルゴリズムが提案されている。以降では、シミュレーテッド分岐アルゴリズムを使って組合せ最適化問題を解く情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラムについて説明する。
　はじめに、シミュレーテッド分岐アルゴリズムの概要について述べる。
　シミュレーテッド分岐アルゴリズムでは、それぞれＮ個ある２つの変数ｘｉ，ｙｉ（ｉ＝１、２、・・・、Ｎ）について、下記の（２）の連立常微分方程式を数値的に解く。Ｎ個の変数ｘｉのそれぞれは、イジングモデルのスピンｓｉに対応している。一方、Ｎ個の変数ｙｉのそれぞれは、運動量に相当している。変数ｘｉ，ｙｉは、いずれも連続変数であるものとする。以下では、変数ｘｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルを第１ベクトル、変数ｙｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルを第２ベクトルとそれぞれよぶものとする。
Figure JPOXMLDOC01-appb-M000002
　ここで、Ｈは、下記の式（３）のハミルトニアンである。
Figure JPOXMLDOC01-appb-M000003
　シミュレーテッド分岐アルゴリズムの時間発展を計算することにより、変数ｘｉ，ｙｉ（ｉ＝１、２、・・・、Ｎ）の値を繰り返し更新し、イジングモデルのスピンｓｉ（ｉ＝１、２、・・・、Ｎ）を求めることができる。すなわち、時間発展の計算が行われる場合、第１ベクトルの要素の値と、第２ベクトルの要素の値が繰り返し更新される。ここでは、時間発展の計算が行われると想定し、各係数を説明する。ただし、シミュレーテッド分岐アルゴリズムでは、必ず時間発展の計算を行わなくてもよい。例えば、後述するようにレプリカ交換法（ｒｅｐｌｉｃａ　ｅｘｃｈａｎｇｅ　ｍｅｔｈｏｄ／ｐａｒａｌｌｅｌ　ｔｅｍｐｅｒｉｎｇ）によってシミュレーテッド分岐アルゴリズムの求解を行ってもよい。
　（２）および（３）において、係数Ｄは、離調（ｄｅｔｕｎｉｎｇ）に相当する。係数ｐ（ｔ）は、ポンピング振幅（ｐｕｍｐｉｎｇ　ａｍｐｌｉｔｕｄｅ）に相当する。シミュレーテッド分岐アルゴリズムの時間発展を計算する場合、係数ｐ（ｔ）は、更新回数に応じて値が単調増加する。係数ｐ（ｔ）の初期値は０に設定されていてもよい。係数Ｋは、正のカー係数（Ｋｅｒｒ　ｃｏｅｆｆｉｃｉｅｎｔ）に相当する。係数ｃとして、定数係数を使うことができる。例えば、係数ｃの値を、シミュレーテッド分岐アルゴリズムによる計算の実行前に決めてもよい。例えば、係数ｃをＪ（２）行列の最大固有値の逆数に近い値に設定することができる。例えば、ｃ＝０．５Ｄ√（Ｎ／２ｎ）という値を使うことができる。ここで、ｎは、組合せ最適化問題に係るグラフのエッジ数である。また、ａ（ｔ）は、時間発展の計算時においてｐ（ｔ）とともに増加する係数である。例えば、ａ（ｔ）として、√（ｐ（ｔ）／Ｋ）を使うことができる。なお、（３）および（４）における局所磁場のベクトルｈｉは、省略すること可能である。
　時間発展の計算を行う場合、係数ｐ（ｔ）の値が所定の値を超えた時に、第１ベクトルにおいて、正値である変数ｘｉを＋１、負値である変数ｘｉを－１にそれぞれ変換すると、スピンｓｉを要素とする解ベクトルを得ることができる。この解ベクトルは、イジング問題の解に相当する。なお、第１ベクトルおよび第２ベクトルの更新回数に基づき、上述の変換を実行し、解ベクトルを求めるか否かの判定を行ってもよい。
　シミュレーテッド分岐アルゴリズムの時間発展を計算する場合、シンプレクティック・オイラー法を使い、上述の（２）を離散的な漸化式に変換し、求解を行うことができる。下記の（４）は、漸化式に変換後のシミュレーテッド分岐アルゴリズムの例を示している。
Figure JPOXMLDOC01-appb-M000004
 　ここで、ｔは、時刻であり、Δｔは、時間ステップ（時間刻み幅）である。なお、（４）では、微分方程式との対応関係を示すために、時刻ｔおよび時間ステップΔｔが使われている。ただし、実際にアルゴリズムをソフトウェアまたはハードウェアに実装する際に必ず時刻ｔおよび時間ステップΔｔが明示的なパラメータとして含まれていなくてもよい。例えば、時間ステップΔｔを１とすれば、実装時のアルゴリズムから時間ステップΔｔを除去することが可能である。アルゴリズムを実装する際に、明示的なパラメータとして時間ｔを含めない場合には、（４）において、ｘｉ（ｔ＋Δｔ）をｘｉ（ｔ）の更新後の値であると解釈すればよい。すなわち、上述の（４）における“ｔ”は、更新前の変数の値、“ｔ＋Δｔ”は、更新後の変数の値を示すものとする。
　シミュレーテッド分岐アルゴリズムの時間発展を計算する場合、ｐ（ｔ）の値を初期値（例えば、０）から所定の値まで増加させた後における変数ｘｉの符号に基づき、スピンｓｉの値を求めることができる。例えば、ｘｉ＞０のときｓｇｎ（ｘｉ）＝１、ｘｉ＜０のときｓｇｎ（ｘｉ）＝－１となる符号関数を使うと、ｐ（ｔ）の値が所定の値まで増加したとき、変数ｘｉを符号関数で変換することによってスピンｓｉの値を求めることができる。符号関数として、例えば、ｘｉ≠０のときに、ｓｇｎ（ｘｉ）＝ｘｉ／｜ｘｉ｜、ｘｉ＝０のときにｓｇｎ（ｘｉ）＝１または－１になる関数を使うことができる。組合せ最適化問題の解（例えば、イジングモデルのスピンｓｉ）を求めるタイミングについては、特に問わない。例えば、第１ベクトルおよび第２ベクトルの更新回数または第１係数ｐの値がしきい値より大きくなったときに組合せ最適化問題の解（解ベクトル）を求めてもよい。
　図６のフローチャートは、時間発展によってシミュレーテッド分岐アルゴリズムの解を計算する場合における処理の例を示している。以下では、図６を参照しながら処理を説明する。
　はじめに、計算サーバは、管理サーバ１より問題に対応する行列Ｊｉｊおよびベクトルｈｉを取得する（ステップＳ１０１）。そして、計算サーバは、係数ｐ（ｔ）およびａ（ｔ）を初期化する（ステップＳ１０２）。例えば、ステップＳ１０２で係数ｐおよびａの値を０にすることができるが、係数ｐおよびａの初期値を限定するものではない。次に、計算サーバは、第１変数ｘｉおよび第２変数ｙｉを初期化する（ステップＳ１０３）。ここで、第１変数ｘｉは、第１ベクトルの要素である。また、第２変数ｙｉは、第２ベクトルの要素である。ステップＳ１０３で計算サーバは、例えば、擬似乱数によってｘｉおよびｙｉを初期化してもよい。ただし、ｘｉおよびｙｉの初期化の方法を限定するものではない。また、これとは異なるタイミングに変数の初期化を行ってもよいし、少なくともいずれかの変数を複数回初期化してもよい。
　次に、計算サーバは、第１ベクトルの要素ｘｉに対応する第２ベクトルの要素ｙｉを重み付け加算することによって第１ベクトルを更新する（ステップＳ１０４）。例えば、ステップＳ１０４では、変数ｘｉにΔｔ×Ｄ×ｙｉを加算することができる。そして、計算サーバは、第２ベクトルの要素ｙｉを更新する（ステップＳ１０５およびＳ１０６）。例えば、ステップＳ１０５では、変数ｙｉにΔｔ×［（ｐ－Ｄ－Ｋ×ｘｉ×ｘｉ）×ｘｉ］を加算することができる。ステップＳ１０６では、さらに変数ｙｉに－Δｔ×ｃ×ｈｉ×ａ－Δｔ×ｃ×ΣＪｉｊ×ｘｊを加算することができる。
　次に、計算サーバは、係数ｐおよびａの値を更新する（ステップＳ１０７）。例えば、係数ｐに一定の値（Δｐ）を加算し、係数ａを更新後の係数ｐの正の平方根に設定することができる。ただし、後述するように、これは係数ｐおよびａの値の更新方法の一例にしかすぎない。そして、計算サーバは、第１ベクトルおよび第２ベクトルの更新回数がしきい値未満であるか否かを判定する（ステップＳ１０８）。更新回数がしきい値未満である場合（ステップＳ１０８のＹＥＳ）、計算サーバは、ステップＳ１０４～Ｓ１０７の処理を再度実行する。更新回数がしきい値以上である場合（ステップＳ１０８のＮＯ）、第１ベクトルの要素ｘｉに基づいて解ベクトルの要素であるスピンｓｉを求める（ステップＳ１０９）。ステップＳ１０９では、例えば、第１ベクトルにおいて、正値である変数ｘｉを＋１、負値である変数ｘｉを－１にそれぞれ変換し、解ベクトルを得ることができる。
　なお、ステップＳ１０８の判定において、更新回数がしきい値未満である場合（ステップＳ１０８のＹＥＳ）に第１ベクトルに基づきハミルトニアンの値を計算し、第１ベクトルおよびハミルトニアンの値を記憶してもよい。これにより、ユーザは、複数の第１ベクトルより最適解に最も近い近似解を選択することが可能となる。
　なお、図６のフローチャートに示した少なくともいずれかの処理を並列的に実行してもよい。例えば、第１ベクトルおよび第２ベクトルのそれぞれが有するＮ個の要素の少なくとも一部が並列的に更新されるよう、ステップＳ１０４～Ｓ１０６の処理を並列的に実行してもよい。例えば、複数台の計算サーバを使って処理を並列化してもよい。複数のプロセッサによって処理を並列化してもよい。ただし、処理の並列化を実現するための実装および処理の並列化の態様を限定するものではない。
　上述のステップＳ１０５～Ｓ１０６に示した変数ｘｉおよびｙｉの更新処理の実行順序は、一例にしかすぎない。したがって、これとは異なる順序で変数ｘｉおよびｙｉの更新処理を実行してもよい。例えば、変数ｘｉの更新処理と変数ｙｉの更新処理が実行される順序が入れ替わっていてもよい。また、各変数の更新処理に含まれるサブ処理の順序も限定しない。例えば、変数ｙｉの更新処理に含まれる加算処理の実行順序が図６の例とは異なっていてもよい。各変数の更新処理を実行するための前提となる処理の実行順序およびタイミングも特に限定しない。例えば、問題項の計算処理が、変数ｘｉの更新処理を含むその他の処理と並行で実行されていてもよい。変数ｘｉおよびｙｉの更新処理、各変数の更新処理に含まれるサブ処理および問題項の計算処理が実行される順序およびタイミングが限定されない点は、以降に示す各フローチャートの処理についても、同様である。
　図７のグラフは、時間発展の計算を行う場合における係数ｐの更新パターンの例を示している。図７の縦軸は、係数ｐの値を示している。また、図７の横軸は、時刻または更新回数に相当する。図７は、初期値ｐ１よりＭ回の更新後の値であるｐＭに至る複数の更新パターンがあることを示している。ｐ１＜ｐＭの条件を満たすのであれば、ｐ１およびｐＭを任意の実数に設定することが可能である。このように、係数ｐの更新パターンは、無数に存在するといえる。
　なお、シミュレーテッド分岐アルゴリズムにおける係数ｐは、シミュレーテッドアニーリング（ＳＡ）における温度Ｔの逆数に相当する。したがって、係数ｐに応じて粒子の遷移確率を調整することができる。係数ｐが小さいほど、粒子の遷移確率が高くなる。一方、係数ｐが大きいほど、粒子の遷移確率が低くなる。このため、シミュレーテッド分岐アルゴリズムにおける係数ｐの更新パターンは、シミュレーテッドアニーリングにおける冷却スケジュールに相当するといえる。
　例えば、図７に示されているように、更新回数に応じて係数ｐを一定の割合でインクリメントしてもよい。また、更新タイミングによって係数ｐを異なる増分でインクリメントしてもよい。さらに、一部の更新タイミングにおいて係数ｐの値を更新し、残りの更新タイミングにおいて係数ｐの値の更新をスキップしてもよい。
　シミュレーテッド分岐アルゴリズムの求解処理においては、複数の局所解を含む解空間のできる限り広い領域を探索し、問題の最適解またはそれに近い近似解に到達することが望ましい。係数ｐおよびａの更新パターンは、探索対象とする解空間の領域と探索の粒度を指定する。このため、大域的な最適解が見つかりやすくなるよう、係数ｐおよびａの更新パターンを調整する必要がある。ただし、一般に求解の対象となるそれぞれの問題ごとに最適な更新パターンを見つけるのは容易であるとは限らない。ユーザが様々な更新パターンを試し、経験的に使用する更新パターンを選ぶことになるケースも想定される。
　そこで、時間発展の計算に代わり、レプリカ交換法によってシミュレーテッド分岐アルゴリズムの求解を行ってもよい。
　図８は、レプリカ交換法によって計算を行う場合に使われる係数ｐを概念的に示したものである。図８の縦軸は、係数ｐの値を示している。レプリカ交換法によって計算を行う場合、情報処理装置は、必ず係数ｐの時間発展を計算しなくてもよい。このため、図８の横軸は、図７との対比のために用意された仮想的な軸となっている。
　図８に示されているように、レプリカ交換法を使う場合、決められた複数の係数ｐの組合せ（ｐ（１）、ｐ（２）、・・・、ｐ（Ｍ））ごとに、第１ベクトルおよび第２ベクトルの更新を行う。以下の説明では、それぞれの係数ｐごとに用意される第１ベクトルおよび第２ベクトルをレプリカとよぶものとする。そして、レプリカ間で第１ベクトルおよび第２ベクトルの交換を行うことにより、解空間内の異なる領域における解の探索を実現することができる。レプリカ交換法では、第１ベクトルおよび第２ベクトルの交換を繰り返すことにより、問題の最適解またはそれに近い近似解に到達する確率を高めることが可能である。
　例えば、計算サーバのプロセッサごとにレプリカを割り当ててもよい。また、情報処理システム内の計算サーバごとにレプリカを割り当ててもよい。また、レプリカによって割り当てられる計算資源の量が異なっていてもよい。レプリカによって割り当てられる計算資源の種類または組合せが異なっていてもよい。ただし、それぞれのレプリカへの計算資源の割り当て方法を限定するものではない。
　レプリカ交換法を使う場合、それぞれのレプリカにおける第１ベクトル、第２ベクトルおよびそれらの関数の計算を並列的に行うことができる（第１レベルの並列化）。また、それぞれが複数の要素を含む第１ベクトルおよび第２ベクトルの更新処理も並列的に行うことができる（第２レベルの並列化）。計算資源において第１レベルの並列化が行われる粒度と、第２レベルの並列化が行われる粒度は異なっていてもよい。前者の粒度より後者の粒度を大きくしてもよい。例えば、第１レベルの並列化を計算サーバ、仮想計算機の粒度で実現することができる。また、第２レベルの並列化をプロセッサの粒度で実現することができる。レプリカ交換法を使うと、複数レベルにおける処理の並列化を行うことができるため、効率的に計算を行うことができる。このため、比較的短時間で問題の最適解またはそれに近い近似解を求めることが可能となる。
［係数ｐの組合せを自動的に決める場合］
　はじめに、レプリカ交換法による計算を行う際に、使用する係数ｐの組合せを自動的に決める場合における処理を説明する。
　例えば、下記の式（５）のハミルトニアンが使われる場合を想定する。
Figure JPOXMLDOC01-appb-M000005
 　ここで、ａ（ｔ）として、例えば、√（ｐ（ｔ））を使うことができる。
　式（５）において、ａ（ｔ）Ｘｉ＝ｘｉおよびａ（ｔ）Ｙｉ＝ｙｉとおき、スケーリングを行うと、下記の式（６）が得られる。
Figure JPOXMLDOC01-appb-M000006
　レプリカ交換法を使う場合、レプリカごとにハミルトニアンの計算が行われる。このため、以下の説明では、対応するレプリカｍ（ｍ＝１、２、・・・、Ｍ）を示すため、下記の式（７）のような表記を用いるものとする。
Figure JPOXMLDOC01-appb-M000007
　統計力学的な平衡分布への収束を保証するために、レプリカｍ１よりレプリカｍ１＋１に遷移する確率Ｗ（ｍ１，ｍ１＋１）と、レプリカｍ１＋１よりレプリカｍ１に遷移する確率Ｗ（ｍ１＋１，ｍ１）は、下記の式（８）の詳細釣り合い条件が課せられる。係数βは、一般的なシミュレーテッドアニーリングにおける逆温度に相当する。レプリカ交換法の各レプリカに係数ｐの組合せ（ｐ（１）、ｐ（２）、・・・、ｐ（Ｍ））を割り当てているため、係数βとして、有限の正定数を使うことができる。
Figure JPOXMLDOC01-appb-M000008
　また、レプリカｍ１とレプリカｍ１＋１との間の交換確率は、下記の式（９）のように表される。
Figure JPOXMLDOC01-appb-M000009
 例えば、データ交換回路５１は、式（９）の値を０以上１未満の範囲の擬似乱数と比較し、前者の値と後者の値との大小関係に基づく条件判定を行うことができる。前者の値より後者の値の方が大きい場合、条件判定は肯定的となり、データ交換回路５１は、レプリカｍ１とレプリカｍ１＋１との間で第１ベクトルおよび第２ベクトルを交換することができる。一方、前者の値より後者の値の方が大きい場合、条件判定は否定的となり、レプリカｍ１とレプリカｍ１＋１との間で第１ベクトルおよび第２ベクトルの交換処理をスキップすることができる。
　例えば、それぞれのレプリカは、演算回路ごとに割り当てられるものであってもよい。この場合、データ交換回路は、演算回路のペアごとに、条件判定を実行し、条件判定が肯定的であるときに演算回路間で計算に使われる第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するように構成されていてもよい。また、データ交換回路は、ペアに含まれる演算回路における第１ベクトル（複数の第１変数）、第２ベクトル（複数の第２変数）および第１係数に基づいて条件判定を実行するように構成されていてもよい。例えば、式（９）を使って条件判定を行ってもよいし、その他の式を使って条件判定を行ってもよい。
　それぞれのレプリカは、計算機ごとに割り当てられるものであってもよい。この場合、情報処理システムの第２計算機は、第１計算機のペアごとに、条件判定を実行し、条件判定が肯定的であるときに第１計算機間で計算に使われる第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するように構成されていてもよい。
　レプリカ交換法では、効率的な計算（平衡状態への収束）を実現するために、係数ｐの値の相違に関わらず交換確率が一定となるよう、それぞれのレプリカにおける係数ｐ（ｍ）の値を設定することが好ましい。特に、複数のレプリカにおける係数ｐ（ｍ）（ｍ＝１、２、・・・、Ｍ）の値が等比数列を形成するよう、各係数ｐ（ｍ）の値を設定すると、レプリカ間における交換確率の値が一定となることが知られている。したがって、それぞれのレプリカにおける係数ｐ（ｍ）として、等比数列を形成する係数ｐ（ｍ）の組合せ（ｍ＝１、２、・・・、Ｍ）を使って計算を行ってもよい。すなわち、それぞれの演算回路で使われる第１係数（例えば、係数ｐ）は、互いに等比数列の関係となっていてもよい。ただし、レプリカ交換法による計算において使われる係数ｐ（ｍ）（ｍ＝１、２、・・・、Ｍ）の値の決め方を限定するものではない。
　各レプリカでは、それぞれＮ個ある２つの変数Ｘ（ｍ） ｉ，Ｙ（ｍ） ｉ（ｉ＝１、２、・・・、Ｎ）について、下記の（１０）の連立常微分方程式を数値的に解く処理が実行される。Ｎ個の変数Ｘ（ｍ） ｉのそれぞれは、イジングモデルのスピンｓｉに対応している。一方、Ｎ個の変数Ｙ（ｍ） ｉのそれぞれは、運動量に相当している。変数Ｘ（ｍ） ｉ，Ｙ（ｍ） ｉは、いずれも連続変数であるものとする。以下では、変数Ｘ（ｍ） ｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルを第１ベクトル、変数Ｙ（ｍ） ｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルを第２ベクトルとそれぞれよぶものとする。
Figure JPOXMLDOC01-appb-M000010
　シンプレクティック・オイラー法を使い、上述の（１０）を離散的な漸化式に変換し、レプリカ交換法によってシミュレーテッド分岐アルゴリズムの計算を行うことができる。下記の（１１）は、漸化式に変換後のシミュレーテッド分岐アルゴリズムの例を示している。
Figure JPOXMLDOC01-appb-M000011
　なお、（１１）のうち、下記の（１２）の項は、イジングエネルギーに由来する。この項の形式は、解きたい問題に応じて決まるため、問題項（ｐｒｏｂｌｅｍ　ｔｅｒｍ）とよぶものとする。
Figure JPOXMLDOC01-appb-M000012
 後述するように、計算において（１２）とは異なる問題項を使ってもよい。
　上述のように、レプリカ交換法によってシミュレーテッド分岐アルゴリズムの計算を行うと、処理の並列度を高め、比較的短時間で問題の最適解またはそれに近い近似解を得ることができる。特に、係数ｐの組合せが自動的に決まる上述の方法を用いると、シミュレーテッド分岐アルゴリズムの実装および実行が容易となる。
　図９のフローチャートは、レプリカ交換法によってシミュレーテッド分岐アルゴリズムの計算を行ったときにおける処理の例を示している。以下では、図９を参照しながら処理を説明する。
　はじめに、計算サーバは、問題に対応する行列Ｊｉｊおよびベクトルｈｉを取得し、カウンタｔおよび係数ｐ（ｍ）、ａ（ｍ）を初期化する（ステップＳ１１１）。例えば、ステップＳ１１１では、等比数列を形成するように、係数ｐ（ｍ）の値を設定することができる。カウンタｔは、例えば、０に設定される。そして、計算サーバは、第１ベクトルの要素Ｘ（ｍ） ｉおよび第２ベクトルの要素Ｙ（ｍ） ｉを初期化し、ｍ１に１を代入する（ステップＳ１１２）。ここで、ｍ１は、レプリカｍの指定に使われるパラメータである。例えば、ステップＳ１１２では、擬似乱数によって第１ベクトルの要素Ｘ（ｍ） ｉおよび第２ベクトルの要素Ｙ（ｍ） ｉの初期値を設定することができる。次に、計算サーバは、第１ベクトルの要素Ｘ（ｍ１） ｉを対応する第２ベクトルの要素Ｙ（ｍ） ｉに基づき更新する（ステップＳ１１３）。例えば、ステップＳ１１３では、第１ベクトルの要素Ｘ（ｍ１） ｉにΔｔ×Ｄ×Ｙ（ｍ１） ｉを加算することができる。そして、計算サーバは、第２ベクトルの要素Ｙ（ｍ１） ｉを更新する（ステップＳ１１４）。例えば、ステップＳ１１４では、第２ベクトルの要素Ｙ（ｍ１） ｉにΔｔ×［（ｐ－Ｄ－Ｋ×Ｘ（ｍ１） ｉ×Ｘ（ｍ１） ｉ）×Ｘ（ｍ１） ｉ］を加算することができる。ステップＳ１１４では、さらに変数Ｙ（ｍ１） ｉに－Δｔ×ｃ×ｈｉ×ａ－Δｔ×ｃ×ΣＪｉｊ×Ｘ（ｍ１） ｉを加算することができる。
　次に、計算サーバは、ｍ１をインクリメントし、ｍ１＜＝Ｍを満たすか否かを判定する（ステップＳ１１５）。ステップＳ１１５の判定が肯定的である場合（ステップＳ１１５のＹＥＳ）、ステップＳ１１３～Ｓ１１５の処理を再度実行する。ステップＳ１１５の判定が否定的である場合（ステップＳ１１５のＮＯ）、計算サーバは、ｋ回ごとに０～Ｍのランダムなｍ１を選択し、遷移条件を満たすのであれば、Ｘ（ｍ１） ｉとＸ（ｍ１＋１） ｉの値、Ｙ（ｍ１） ｉとＹ（ｍ１＋１） ｉの値をそれぞれ交換する（ステップＳ１１６）。例えば、計算サーバは、０以上１未満の擬似乱数ｒを生成し、ｒの値を上述の式（９）の値と比較することができる。上述の式（９）の値の方が擬似乱数ｒの値より大きい場合、遷移条件を満たしていると判定することができる。一方、擬似乱数ｒの値が上述の式（９）の値以上である場合、遷移条件が満たされていないと判定することができる。そして、計算サーバは、カウンタｔを更新する（ステップＳ１１７）。例えば、ステップＳ１１７では、ｔにΔｔを加算することができる。
　次に、カウンタｔがＴ未満であるか否かを判定する（ステップＳ１１８）。ステップＳ１１８の判定が肯定的である場合（ステップＳ１１８のＹＥＳ）、計算サーバは、ステップＳ１１３～Ｓ１１７の処理を再度実行する。ステップＳ１１８の判定が否定的である場合（ステップＳ１１８のＮＯ）、計算サーバは、第１ベクトルの要素（変数）Ｘ（ｍ） ｉの値に基づいて解ベクトルの要素であるスピンｓｉを求める（ステップＳ１１９）。例えば、ステップＳ１１９では、正値である変数Ｘ（ｍ） ｉを＋１、負値である変数Ｘ（ｍ） ｉを－１にそれぞれ変換すればよい。なお、ステップＳ１１８では、代わりに、ステップＳ１１６の処理の実行回数がしきい値を超えたか否かの判定を行ってもよい。
　ステップＳ１１９における変数Ｘ（ｍ） ｉとして、１～Ｍのいずれかのレプリカより選択された変数Ｘ（ｍ１） ｉの代表値を使うことができる。ここで、代表値は、１～Ｍのいずれかのレプリカより選択された変数Ｘ（ｍ１） ｉであってもよいし、１～Ｍの少なくともいずれかのレプリカより選択された変数Ｘ（ｍ１） ｉの平均値であってもよい。
　上述の図９のフローチャートでは、ステップＳ１１２～Ｓ１１４およびＳ１１６のすくなくとも一部において第１ベクトルの要素および第２ベクトルの要素（ｉ＝１，２，・・・Ｎ）の更新処理を並列的に実行してもよい（第２レベルの並列化）。また、それぞれのレプリカに係る処理の少なくとも一部を並列化してもよい（第１レベルの並列化）。例えば、計算サーバ、仮想計算機またはプロセッサごとに異なるレプリカに係るステップＳ１１２～Ｓ１１４の処理を実行させてもよい。この場合、ステップＳ１１６の処理は、いずれかの計算サーバ上で動作するプログラムによって実行されてもよいし、上述の管理サーバ１の制御部１３によって実行されてもよい。また、ステップＳ１１９の処理は、管理サーバ１など、その他のサーバによって実行されてもよい。
　ここでは、本実施形態による情報処理方法、記憶媒体およびプログラムについて述べる。
　情報処理方法は、複数の第１計算機と、第２計算機とを使うものであってもよい。例えば、複数の第１計算機は、それぞれが第１変数を要素とする第１ベクトルおよび第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている。例えば、第２計算機は、第１計算機間におけるデータの交換を実行するように構成されている。この場合、情報処理方法は、それぞれの第１計算機が第１変数を対応する第２変数に基づいて更新するステップと、それぞれの第１計算機が第１変数を第１係数で重み付けし対応する第２変数に加算するステップと、それぞれの第１計算機が複数の第１変数を用いて問題項を計算するステップと、それぞれの第１計算機が問題項を第２変数に加算するステップと、第２計算機が第１計算機間で第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するステップとを含んでいてもよい。
　また、情報処理方法は、複数の第１プロセッサと、第２プロセッサとを使うものであってもよい。例えば、複数の第１プロセッサは、それぞれが第１変数を要素とする第１ベクトルおよび第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている。例えば、第２プロセッサは、第１プロセッサ間におけるデータの交換を実行するように構成されている。この場合、情報処理方法は、それぞれの第１プロセッサが第１変数を対応する第２変数に基づいて更新するステップと、それぞれの第１プロセッサが第１変数を第１係数で重み付けし対応する第２変数に加算するステップと、それぞれの第１プロセッサが複数の第１変数を用いて問題項を計算するステップと、それぞれの第１プロセッサが問題項を第２変数に加算するステップと、第２プロセッサが第１プロセッサ間で第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するステップとを含んでいてもよい。本実施形態によるプログラムは、当該情報処理方法と同様のステップを含んでいてもよい。さらに、本実施形態による非一時的なコンピュータ可読な記憶媒体は、当該プログラムを格納するものであってもよい。
［係数ｐの値が変動する場合］
　ここまでは、各レプリカにおける係数ｐ（ｍ）、ａ（ｍ）の値がループ処理の前のステップ（図９のステップＳ１１１）で決まり、以降のループ処理を含むステップでは、固定された係数ｐ（ｍ）、ａ（ｍ）の値が使われる場合を想定し、説明を行った。ただし、レプリカ交換法によってシミュレーテッド分岐アルゴリズムの計算が行われる場合、必ず係数ｐ（ｍ）、ａ（ｍ）の値は、固定されていなくてもよい。
　例えば、情報処理装置は、下記の式（１３）に基づいて係数ｐ（ｍ）の値を変動させてもよい。式（１３）は、レプリカの交換タイミングまたは交換判定タイミングによって変動しうる係数ｐ（ｍ）の一例である。
Figure JPOXMLDOC01-appb-M000013
 式（１３）のｐ（ｍ）（ｔ）は、係数ｐ（ｍ）の値がパラメータｔに依存していることを示している。
条件（ａ）
はじめに、パラメータｔの値が０またはＴ０の偶数倍である場合におけるｐ（ｍ）の値について説明する。
　レプリカ番号ｍが偶数である場合、式（１３）のｃｏｓ（・・・）の値は、＋１となり、ｐ（ｍ）＝ｐ０ （ｍ＋１）が成り立つ。一方、レプリカ番号ｍが奇数である場合、式（１３）のｃｏｓ（・・・）の値は、－１となり、ｐ（ｍ）＝ｐ０ （ｍ）が成り立つ。
条件（ｂ）
　次に、パラメータｔの値がＴ０の奇数倍である場合におけるｐ（ｍ）の値について説明する。
　レプリカ番号ｍが偶数である場合、式（１３）のｃｏｓ（・・・）の値は、－１となり、ｐ（ｍ）＝ｐ０ （ｍ）が成り立つ。一方、レプリカ番号ｍが奇数である場合、式（１３）のｃｏｓ（・・・）の値は、＋１となり、ｐ（ｍ）＝ｐ０ （ｍ＋１）が成り立つ。
条件（ｃ）
　最後に、パラメータｔの値について、上述の条件（ａ）および（ｂ）のいずれにも該当しない場合におけるｐ（ｍ）の値について説明する。
　この場合、式（１３）のｃｏｓ（・・・）は、（－１，＋１）の範囲の値をとるため、ｐ（ｍ）は、αｐ０ （ｍ＋１）＋βｐ０ （ｍ）のようなｐ０ （ｍ＋１）とｐ０ （ｍ）がそれぞれ重み付けされた値をとる。ｃｏｓ（・・・）の周期性より、番号ｍが偶数のレプリカと、番号ｍが奇数のレプリカで異なる値のｐ（ｍ）が使われることになる。このように、それぞれの演算回路は、第１ベクトルおよび第２ベクトルの更新回数に応じて、第１係数の値を周期的に変化させてもよい。
　上述の式（１３）で定義される係数ｐ（ｍ）を使うことにより、レプリカ間の交換判定の行われるタイミングにおいて、下記の式（１４）が成り立つ。
Figure JPOXMLDOC01-appb-M000014
 レプリカ間の交換判定タイミングにおいて、式（１４）の関係が成立する場合、それぞれの交換判定タイミングにおいてレプリカｍ１とレプリカｍ１＋１について、第１ベクトルおよび第２ベクトルの交換が行われる。それぞれの交換判定タイミングにおいて式（９）に示した計算を行うことが不要となるため、計算資源の消費を抑制することができる。また、図９のフローチャートにおける交換判定タイミング（ステップＳ１１６）で、交換を行うか否かの判定処理をスキップすることができるため、図９の交換判定タイミングは、交換タイミングに置き換わる。
　式（１３）に示した係数ｐ（ｍ）を使う場合、ユーザは、事前にｐ０ （ｍ＋１）およびｐ０ （ｍ）の値を決める必要がある。しかし、時間発展の計算を行う場合のように係数ｐ（ｍ）の更新パターンおよび係数ａ（ｍ）の定義を決める必要がないため、シミュレーテッド分岐アルゴリズムの実行に要する工数を削減することが可能である。また、式（１３）に示した係数ｐ（ｍ）を使うことにより、レプリカ間における第１ベクトルおよび第２ベクトルの交換を効率的に行い、比較的短時間で問題の最適解またはそれに近い近似解を得ることが可能となる。
　なお、式（１３）に示した係数ｐ（ｍ）の定義は、一例にしかすぎない。したがって、これとは異なる定義の係数ｐ（ｍ）を用いることを妨げるものではない。
　図１０のフローチャートは、変動する係数ｐ（ｍ）を使って計算を行ったときにおける処理の例を示している。以下では、図１０を参照しながら処理を説明する。
　はじめに、計算サーバは、問題に対応する行列Ｊｉｊおよびベクトルｈｉを取得し、カウンタｔおよび係数ｐ（ｍ）、ａ（ｍ）を初期化する（ステップＳ１２１）。カウンタｔは、例えば、０に設定される。式（１３）の定義にしたがう場合、番号ｍが偶数のレプリカでは、ｐ（ｍ）＝ｐ０ （ｍ＋１）に設定される。また、番号ｍが奇数のレプリカでは、ｐ（ｍ）＝ｐ０ （ｍ）に設定される。
　ステップＳ１２２～ステップＳ１２５の処理は、それぞれ図９のステップＳ１１２～ステップＳ１１５の処理に相当するため、詳細な説明を省略する。ステップＳ１２５の判定が肯定的である場合（ステップＳ１２５のＹＥＳ）、計算サーバは、ｋ回ごとに０～Ｍのランダムなｍ１を選択し、Ｘ（ｍ１） ｉとＸ（ｍ１＋１） ｉの値、Ｙ（ｍ１） ｉとＹ（ｍ１＋１） ｉの値をそれぞれ交換する（ステップＳ１２６）。式（１４）の関係が成立する場合には、レプリカ間で交換を行うか否かの判定を行う必要がなくなるため、ステップＳ１２６は、レプリカ間における第１ベクトルおよび第２ベクトルの交換タイミングに相当する。そして、計算サーバは、カウンタｔを更新する（ステップＳ１２７）。例えば、ステップＳ１２７では、ｔにΔｔを加算することができる。
　次に、カウンタｔがＴ未満であるか否かを判定する（ステップＳ１２８）。ステップＳ１２８の判定が肯定的である場合（ステップＳ１２８のＹＥＳ）、計算サーバは、ステップＳ１２３～Ｓ１２７の処理を再度実行する。ステップＳ１２８の判定が否定的である場合（ステップＳ１２８のＮＯ）、計算サーバは、第１ベクトルの要素（変数）Ｘ（ｍ） ｉの値に基づいて解ベクトルの要素であるスピンｓｉを求める（ステップＳ１２９）。例えば、ステップＳ１２９では、正値である変数Ｘ（ｍ） ｉを＋１、負値である変数Ｘ（ｍ） ｉを－１にそれぞれ変換すればよい。ステップＳ１２９における変数Ｘ（ｍ） ｉとして、使うことができる代表値の選び方および計算方法を限定しない点は、図９のステップＳ１１９と同様である。また、ステップＳ１２９の処理は、計算サーバ以外のサーバ（計算機）によって実行されてもよい。
　上述の図１０のフローチャートでは、ステップＳ１２２～Ｓ１２４およびＳ１２６のすくなくとも一部において第１ベクトルの要素および第２ベクトルの要素（ｉ＝１，２，・・・Ｎ）の更新処理を並列的に実行してもよい（第２レベルの並列化）。また、それぞれのレプリカに係る処理の少なくとも一部を並列化してもよい（第１レベルの並列化）。例えば、計算サーバ、仮想計算機またはプロセッサごとに異なるレプリカに係るステップＳ１２２～Ｓ１２４の処理を実行させてもよい。この場合、ステップＳ１２６の処理は、いずれかの計算サーバ上で動作するプログラムによって実行されてもよいし、上述の管理サーバ１の制御部１３によって実行されてもよい。
　情報処理装置が式（１３）に示した係数ｐ（ｍ）を使い、図１０のフローチャートの処理を実行することにより、計算量を抑えつつ、処理の並列度を高め、比較的短時間で問題の最適解またはそれに近い近似解を得ることが可能となる。
［レプリカ交換処理の並列化］
　レプリカ交換法によるシミュレーテッド分岐アルゴリズムの計算において、レプリカ間における第１ベクトルおよび第２ベクトルの交換処理を行うためには、少なくともレプリカ交換タイミングにおいて、交換対象のレプリカｍ１およびｍ１＋１の処理が同期している必要がある。このため、実装次第でレプリカ交換処理は、レプリカ交換法によるシミュレーテッド分岐アルゴリズムの高速化におけるボトルネックとなりうる。そこで、レプリカ交換処理を並列に実行し、交換対象のレプリカ間で処理の同期をとりやすくする対策がとることができる。
　例えば、偶数番号のレプリカｍ１と、レプリカｍ１＋１との間における第１ベクトルおよび第２ベクトルの交換処理（第１交換処理）と、奇数番号のレプリカｍ１と、レプリカｍ１＋１との間における第１ベクトルおよび第２ベクトルの交換処理（第２交換処理）とを交互に実行することができる。これにより、レプリカ間における第１ベクトルおよび第２ベクトルの交換処理を並列化することが可能となる。
　図１１のフローチャートは、並列的にレプリカ交換を行う場合における処理の例を示している。以下では、図１１を参照しながら処理を説明する。
　はじめに、計算サーバは、問題に対応する行列Ｊｉｊおよびベクトルｈｉを取得し、カウンタｔおよび係数ｐ（ｍ）、ａ（ｍ）を初期化する（ステップＳ１３１）。例えば、ステップＳ１３１における係数ｐ（ｍ）の初期化の方法は、上述のステップＳ１２１と同様であってもよい。また、ステップＳ１３１における係数ｐ（ｍ）の初期化の方法は、上述のステップＳ１１１と同様であってもよい。また、ステップＳ１３１では、経験側で得られた係数ｐ（ｍ）の組合せを使ってもよい。カウンタｔは、例えば、０に設定される。
　そして、計算サーバは、第１ベクトルの要素Ｘ（ｍ） ｉおよび第２ベクトルの要素Ｙ（ｍ） ｉを初期化し、ｍ１に１を代入する（ステップＳ１３２）。ここで、ｍ１は、レプリカｍの指定に使われるパラメータである。例えば、ステップＳ１３２では、擬似乱数によって第１ベクトルの要素Ｘ（ｍ） ｉおよび第２ベクトルの要素Ｙ（ｍ） ｉの初期値を設定することができる。次に、計算サーバは、第１ベクトルの要素Ｘ（ｍ１） ｉを対応する第２ベクトルの要素Ｙ（ｍ） ｉに基づき更新する（ステップＳ１３３）。例えば、ステップＳ１３３では、第１ベクトルの要素Ｘ（ｍ１） ｉにΔｔ×Ｄ×Ｙ（ｍ１） ｉを加算することができる。そして、計算サーバは、第２ベクトルの要素Ｙ（ｍ１） ｉを更新する（ステップＳ１３４）。例えば、ステップＳ１３４では、第２ベクトルの要素Ｙ（ｍ１） ｉにΔｔ×［（ｐ－Ｄ－Ｋ×Ｘ（ｍ１） ｉ×Ｘ（ｍ１） ｉ）×Ｘ（ｍ１） ｉ］を加算することができる。ステップＳ１３４では、さらに第２ベクトルの要素Ｙ（ｍ１） ｉに－Δｔ×ｃ×ｈｉ×ａ－Δｔ×ｃ×ΣＪｉｊ×Ｘ（ｍ１） ｉを加算することができる。
　次に、計算サーバは、ｋ回ごとに偶数であるｍ１または奇数であるｍ１を交互に選択して遷移条件を判定し、遷移条件が満たされるのであれば、Ｘ（ｍ１） ｉとＸ（ｍ１＋１） ｉの値、Ｙ（ｍ１） ｉとＹ（ｍ１＋１） ｉの値をそれぞれ交換する（ステップＳ１３５）。例えば、計算サーバは、０以上１未満の擬似乱数ｒを生成し、ｒの値を上述の式（９）の値と比較することができる。上述の式（９）の値の方が擬似乱数ｒの値より大きい場合、遷移条件を満たしていると判定することができる。一方、擬似乱数ｒの値が上述の式（９）の値以上である場合、遷移条件が満たされていないと判定することができる。
　そして、計算サーバは、ｍ１をインクリメントし、ｍ１＜＝Ｍを満たすか否かを判定する（ステップＳ１３６）。ステップＳ１３６の判定が肯定的である場合（ステップＳ１３６のＹＥＳ）、ステップＳ１３３～Ｓ１３５の処理を再度実行する。ステップＳ１３６の判定が否定的である場合（ステップＳ１３６のＮＯ）、計算サーバは、カウンタｔを更新する（ステップＳ１３７）。例えば、ステップＳ１３７では、ｔにΔｔを加算することができる。
　次に、カウンタｔがＴ未満であるか否かを判定する（ステップＳ１３８）。ステップＳ１３８の判定が肯定的である場合（ステップＳ１３８のＹＥＳ）、計算サーバは、ステップＳ１３３～Ｓ１３７の処理を再度実行する。ステップＳ１３８の判定が否定的である場合（ステップＳ１３８のＮＯ）、計算サーバは、第１ベクトルの要素（変数）Ｘ（ｍ） ｉの値に基づいて解ベクトルの要素であるスピンｓｉを求める（ステップＳ１３９）。例えば、ステップＳ１３９では、正値である変数Ｘ（ｍ） ｉを＋１、負値である変数Ｘ（ｍ） ｉを－１にそれぞれ変換すればよい。ステップＳ１３９における変数Ｘ（ｍ） ｉとして、使うことができる代表値の選び方および計算方法を限定しない点は、図９のステップＳ１１９と同様である。また、ステップＳ１３９の処理は、計算サーバ以外のサーバ（情報処理装置）によって実行されてもよい。
　上述の図１１のフローチャートでは、ステップＳ１３２～Ｓ１３５のすくなくとも一部において第１ベクトルの要素および第２ベクトルの要素（ｉ＝１，２，・・・Ｎ）の更新処理を並列的に実行してもよい（第２レベルの並列化）。また、それぞれのレプリカに係る処理の少なくとも一部を並列化してもよい（第１レベルの並列化）。例えば、計算サーバ、仮想計算機またはプロセッサごとに異なるレプリカに係るステッＳ１３２～Ｓ１３５の処理を実行させてもよい。この場合、ステップＳ１３５の処理は、いずれかの計算サーバ上で動作するプログラムによって実行されてもよいし、上述の管理サーバ１の制御部１３によって実行されてもよい。
　図１２は、レプリカ交換法における処理フローを概念的に示している。図１２において、それぞれの正方形は、番号が１～Ｍのレプリカに相当する。
　図１２のステップＳ１０は、各レプリカについて第１ベクトルおよび第２ベクトルの更新処理が実行されているときを示している。そして、ステップＳ１１では、偶数番号のレプリカについて、遷移条件を満たすか否かを判定し、遷移条件が満たされる場合には、第１ベクトルおよび第２ベクトルの交換処理が実行されている。また、ステップＳ１２では、偶数番号のレプリカについて、遷移条件を満たすか否かを判定し、遷移条件が満たされる場合には、第１ベクトルおよび第２ベクトルの交換処理が実行されている。遷移条件の判定は、上述のように、０以上１未満の擬似乱数ｒを生成し、ｒの値を上述の式（９）の値と比較することによって行うことができる。
　データ交換回路は、タイミングによって異なる演算回路のペアを選択し、第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するように構成されていてもよい。また、データ交換回路は、複数の演算回路において、第１組合せのペアと、第２組合せのペアとを交互に選択し、第１ベクトルおよび第２ベクトルの交換、または第１係数の交換の少なくともいずれかを実行するように構成されていてもよい。例えば、番号１～番号６のレプリカがある場合を想定する。番号１～番号６のレプリカは、それぞれいずれかの演算回路に割り当てられているものとする。この場合、第１組合せには、番号１および番号２のレプリカのペアと、番号３および番号４のレプリカのペアと、番号５および番号６のレプリカのペアとが含まれていてもよい。一方、第２組合せには、番号２および番号３のレプリカのペアと、番号４および番号５のレプリカのペアと、番号６および番号１のレプリカのペアとが含まれていてもよい。
　図１２に示されているように、タイミングによってステップＳ１１の処理と、ステップＳ１２の処理を交互に実施することができる。これにより、レプリカ交換処理を並列的に実行し、一層短時間で問題の最適解またはそれに近い近似解を得ることが可能となる。
　図１３は、レプリカ交換法に関連する式の例を示している。図１３の式（Ｅ１）は、それぞれのレプリカで使われるハミルトニアンの例を示している。式（Ｅ２）は、カノニカルアンサンブルの例を示している。式（Ｅ３）は、一般的な詳細釣り合いの式を示している。式（Ｅ４）は、カノニカルアンサンブルにおける詳細釣り合いの式を示している。式（Ｅ５）は、カノニカルアンサンブルにおける詳細釣り合いの式の項を整理したものである。式（Ｅ６）は、レプリカｍ１とレプリカｍ２の間における遷移確率を示している。
　図１４のグラフは、レプリカ交換法によってシミュレーテッド分岐アルゴリズムの計算が行われる場合と、時間発展によってシミュレーテッド分岐アルゴリズムの計算が行われる場合のそれぞれにおいて、最適解を得られるまでに必要な試行回数の例を示したヒストグラムである。図１４は、シミュレーテッド分岐アルゴリズムを使って４８ノード９６エッジのハミルトン閉路問題を解いた場合におけるデータを示している。図１４の縦軸は、最適解が得られた回数を示している。レプリカ交換法による計算では、１０個のレプリカを使って計算を行っている。したがって、レプリカ交換法による１回の計算における試行回数は、１０回とカウントしている。なお、時間発展による計算の場合には、２００サイクルごとに最適解が得られているか否かの判定を行っている。一方、レプリカ交換法による計算の場合には、５０サイクルごとに最適解が得られているか否かの判定を行っている。図１４の結果を参照すると、レプリカ交換法を使うと、時間発展を計算する場合と比べて少ない試行回数で最適解を得られている。
　レプリカ交換法による計算を行うことにより、情報処理システムに冗長性を持たせることができる。例えば、いずれかのレプリカにおける第１ベクトルおよび第２ベクトルの計算を実行している計算機が異常停止したとしても、その他のレプリカにおける第１ベクトルおよび第２ベクトルの計算を実行している計算機が稼働を続けているのであれば、情報処理システムにおける計算処理を継続させることができる。この場合、その他のレプリカにおけるデータに基づき、異常停止した計算機で使われていたレプリカのデータを復旧させてもよい。いずれかのレプリカにおける第１ベクトルおよび第２ベクトルの計算を実行しているプロセッサが異常停止した場合についても、その他のプロセッサにおいて計算を継続することが可能である。
［多体相互作用の項を含む計算］
　シミュレーテッド分岐アルゴリズムを使うことにより、３次以上の目的関数を有する組合せ最適化問題を解くことも可能である。２値変数を変数とする３次以上の目的関数を最小化する変数の組合せを求める問題は、ＨＯＢＯ（Ｈｉｇｈｅｒ　Ｏｒｄｅｒ　Ｂｉｎａｒｙ　Ｏｐｔｉｍｉｚａｔｉｏｎ）問題とよばれる。ＨＯＢＯ問題を扱う場合、高次へ拡張されたイジングモデルにおけるエネルギー式として、下記の式（１５）を使うことができる。
Figure JPOXMLDOC01-appb-M000015
 ここで、Ｊ（ｎ）はｎ階テンソルであり、式（１）の局所磁場ｈｉと結合係数の行列Ｊを一般化させたものである。例えば、テンソルＪ（１）は、局所磁場ｈｉのベクトルに相当する。ｎ階テンソルＪ（ｎ）では、複数の添え字に同じ値があるとき、要素の値は０となる。式（１５）では、３次の項までが示されているが、それより高次の項も式（１５）と同様に定義することができる。式（１５）は多体相互作用を含むイジングモデルのエネルギーに相当する。
　なお、ＱＵＢＯと、ＨＯＢＯはいずれも、制約なし多項式２値変数最適化（ＰＵＢＯ：Ｐｏｌｙｎｏｍｉａｌ　Ｕｎｃｏｎｓｔｒａｉｎｅｄ　Ｂｉｎａｒｙ　Ｏｐｔｉｍｉｚａｔｉｏｎ）の１種であるといえる。すなわち、ＰＵＢＯのうち、２次の目的関数を有する組合せ最適化問題は、ＱＵＢＯである。また、ＰＵＢＯのうち、３次以上の目的関数を有する組合せ最適化問題は、ＨＯＢＯであるといえる。
　シミュレーテッド分岐アルゴリズムを使ってＨＯＢＯ問題を解く場合、上述の式（５）のハミルトニアンＨを下記の式（１６）のハミルトニアンＨに置き換えればよい。
Figure JPOXMLDOC01-appb-M000016
　また、式（１６）より下記の式（１７）に示した問題項が導かれる。
Figure JPOXMLDOC01-appb-M000017
 （１７）の問題項ｚｉは、（１６）の２番目の式を、いずれかの変数ｘｉ（第１ベクトルの要素）について偏微分した形式をとっている。偏微分される変数ｘｉは、インデックスｉによって異なる。ここで、変数ｘｉのインデックスｉは、第１ベクトルの要素および第２ベクトルの要素を指定するインデックスに相当する。
　多体相互作用の項を含む計算を行う場合、上述の（１１）の漸化式は、下記の（１８）の漸化式に置き換わる。
Figure JPOXMLDOC01-appb-M000018
 　（１８）は、（１１）の漸化式をさらに一般化したものに相当する。
　上述に示した問題項は、本実施形態による情報処理装置が使うことができる問題項の例にしかすぎない。したがって、計算で使われる問題項の形式は、これらとは異なるものであってもよい。
［アルゴリズムの変形例］
　ここでは、シミュレーテッド分岐アルゴリズムの変形例について説明する。例えば、誤差の軽減または計算時間の短縮を目的に、上述のシミュレーテッド分岐アルゴリズムに各種の変形を行ってもよい。
　例えば、計算の誤差を軽減するために、第１ベクトルの要素の更新時に追加の処理を実行してもよい。例えば、更新によって第１ベクトルの要素Ｘ（ｎ） ｉの絶対値が１より大きくなったとき、第１ベクトルの要素Ｘ（ｎ） ｉの値をｓｇｎ（Ｘ（ｎ） ｉ）に置き換える。すなわち、更新によってＸ（ｎ） ｉ＞１となったとき、変数Ｘ（ｎ） ｉの値は１に設定される。また、更新によってＸ（ｎ） ｉ＜－１となったとき、変数Ｘ（ｎ） ｉの値は－１に設定される。これにより、変数Ｘ（ｎ） ｉを使ってスピンｓｉをより高い精度で近似することが可能となる。このような処理を含めることにより、アルゴリズムは、ｘｉ＝±１の位置に壁があるＮ粒子の物理モデルと等価になる。より一般的に述べると、演算回路は、値が第１値より小さい第１変数を第１値に設定し、値が第２値より大きい第１変数を第２値に設定するように構成されていてもよい。
　さらに、更新によってＸ（ｎ） ｉ＞１となったとき、変数Ｘ（ｎ） ｉに対応する変数Ｙ（ｎ） ｉに係数ｒｆを乗算してもよい。例えば、－１＜ｒ≦０の係数ｒｆを使うと、上記の壁は、反射係数ｒｆの壁となる。特に、ｒｆ＝０の係数ｒｆを使った場合、アルゴリズムは、Ｘ（ｎ） ｉ＝±１の位置に完全非弾性衝突の起こる壁がある物理モデルと等価になる。より一般的に述べると、演算回路は、値が第１値より小さい第１変数に対応する第２変数、または、第２値より大きい第１変数に対応する第２変数を、もとの第２変数に、第２係数を乗じた値に更新するように構成されていてもよい。例えば、演算回路は、値が－１より小さい第１変数に対応する第２変数、または、値が１より大きい第１変数に対応する第２変数を、もとの第２変数に第２係数を乗じた値に更新するように構成されていてもよい。ここで、第２係数は上述の係数ｒｆに相当する。
　なお、演算回路は、更新によってＸ（ｎ） ｉ＞１となったとき、変数Ｘ（ｎ） ｉに対応する変数Ｙ（ｎ） ｉの値を擬似乱数に設定してもよい。例えば、［－０．１，０．１］の範囲の乱数を使うことができる。すなわち、演算回路は、値が第１値より小さい第１変数に対応する第２変数の値、または、値が第２値より大きい第１変数に対応する第２変数の値を、擬似乱数に設定するように構成されていてもよい。
　以上のようにして｜Ｘ（ｎ） ｉ｜＞１となることを抑止するように更新処理を実行すれば、（１１）および（１８）の非線形項Ｋ×Ｘ（ｎ１） ｉ（ｔ＋Δｔ）×Ｘ（ｎ１） ｉ（ｔ＋Δｔ）を除去しても、Ｘ（ｎ） ｉの値が発散することはなくなる。したがって、下記の（１９）に示したアルゴリズムを使うことが可能となる。
Figure JPOXMLDOC01-appb-M000019
　（１９）のアルゴリズムでは、問題項において、離散変数ではなく、連続変数Ｘが使われている。このため、本来の組合せ最適化問題で使われている離散変数との誤差が生ずる可能性がある。この誤差を軽減するために、下記の（２０）のように、問題項の計算において、連続変数Ｘの代わりに、連続変数Ｘを符号関数で変換した値ｓｇｎ（Ｘ）を使うことができる。
Figure JPOXMLDOC01-appb-M000020
 　（２０）において、ｓｇｎ（Ｘ）は、スピンｓに相当する。
　（２０）では、問題項の中の１階のテンソルを含む項の係数αを定数（例えば、α＝１）にしてもよい。（２０）のアルゴリズムでは、問題項で現れるスピンどうしの積が必ず－１または１のいずれかの値をとるため、高次の目的関数を有するＨＯＭＯ問題を扱った場合、積演算による誤差の発生を防ぐことができる。上述の（２０）のアルゴリズムのように、計算サーバが計算するデータは、さらに、変数ｓｉ（ｉ＝１、２、・・・、Ｎ）を要素とするスピンのベクトル（ｓ１，ｓ２，・・・，ｓＮ）を含んでいてもよい。第１ベクトルのそれぞれの要素を符号関数で変換することにより、スピンのベクトルを得ることができる。
［変数の更新処理の並列化の例］
　以下では、シミュレーテッド分岐アルゴリズムの計算時における変数の更新処理の並列化の例について説明する。
　はじめに、ＰＣクラスタへシミュレーテッド分岐アルゴリズムを実装した例について説明する。ＰＣクラスタとは、複数台のコンピュータを接続し、１台のコンピュータでは得られない計算性能を実現するシステムである。例えば、図１に示した情報処理システム１００は、複数台の計算サーバおよびプロセッサを含んでおり、ＰＣクラスタとして利用することが可能である。例えば、ＰＣクラスタにおいては、ＭＰＩ（Ｍｅｓｓａｇｅ　Ｐａｓｓｉｎｇ　Ｉｎｔｅｒｆａｃｅ）を使うことにより、情報処理システム１００のような複数の計算サーバにメモリが分散して配置されている構成でも並列的な計算を実行することが可能である。例えば、ＭＰＩを使って管理サーバ１の制御プログラム１４Ｅ、各計算サーバの計算プログラム３４Ｂおよび制御プログラム３４Ｃを実装することができる。
　ＰＣクラスタで利用するプロセッサ数がＱである場合、それぞれのプロセッサに、第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）に含まれる変数ｘｉのうち、Ｌ個の変数の計算を行わせることができる。同様に、それぞれのプロセッサに、第２ベクトル（ｙ１，ｙ２，・・・，ｙＮ）に含まれる変数ｙｉのうち、Ｌ個の変数の計算を行わせることができる。すなわち、プロセッサ＃ｊ（ｊ＝１，２，・・・，Ｑ）は、変数｛ｘｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝および｛ｙｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の計算を行う。また、プロセッサ＃ｊによる｛ｙｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の計算に必要な下記の（２１）に示されたテンソルＪ（ｎ）は、プロセッサ＃ｊがアクセス可能な記憶領域（例えば、レジスタ、キャッシュ、メモリなど）に保存されるものとする。
Figure JPOXMLDOC01-appb-M000021
　ここでは、それぞれのプロセッサが第１ベクトルおよび第２ベクトルの一定数の変数を計算する場合を説明した。ただし、プロセッサによって、計算する第１ベクトルおよび第２ベクトルの変数の数が異なっていてもよい。例えば、計算サーバに実装されるプロセッサによって性能差がある場合、プロセッサの性能に応じて計算対象とする変数の数を決めることができる。
　変数ｙｉの値を更新するためには、第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）のすべての成分の値が必要となる。２値変数への変換は、例えば、符号関数ｓｇｎ（）を使うことによって行うことができる。そこで、Ａｌｌｇａｔｈｅｒ関数を使い、第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）のすべての成分の値をＱ個のプロセッサに共有させることができる。第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）については、プロセッサ間での値の共有が必要であるものの、第２ベクトル（ｙ１，ｙ２，・・・，ｙＮ）およびテンソルＪ（ｎ）については、プロセッサ間での値の共有を行うことは必須ではない。プロセッサ間でのデータの共有は、例えば、プロセッサ間通信を使ったり、共有メモリにデータを保存したりすることによって実現することができる。
　プロセッサ＃ｊは、問題項｛ｚｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の値を計算する。そして、プロセッサ＃ｊは、計算した問題項｛｛ｚｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の値に基づき、変数｛ｙｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝を更新する。
　上述の各式に示したように、問題項のベクトル（ｚ１，ｚ２，・・・，ｚＮ）の計算では、テンソルＪ（ｎ）と、ベクトル（ｘ１，ｘ２，・・・，ｘＮ）との積の計算を含む、積和演算が必要である。積和演算は、上述のアルゴリズムにおいて最も計算量の大きい処理であり、計算速度の向上においてボトルネックとなりうる。そこで、ＰＣクラスタの実装では、積和演算を、Ｑ＝Ｎ／Ｌ個のプロセッサに分散して並列的に実行し、計算時間の短縮をはかることができる。
　図１５は、マルチプロセッサ構成の例を概略的に示している。図１５の複数の計算ノードは、例えば、情報処理システム１００の複数の計算サーバに相当する。また、図１５の高速リンクは、例えば、情報処理システム１００のケーブル４ａ～４ｃおよびスイッチ５によって形成された計算サーバ間のインターコネクトに相当する。図１５の共有メモリは、例えば、共有メモリ３２に相当する。図１５のプロセッサは、例えば、各計算サーバのプロセッサ３３Ａ～３３Ｄに相当している。なお、図１５には複数の計算ノードが示されているが、単一計算ノードの構成を用いることを妨げるものではない。
　図１５には、各構成要素に配置されるデータおよび構成要素間で転送されるデータが示されている。各プロセッサでは、変数ｘｉ、ｙｉの値が計算される。また、プロセッサと共有メモリ間では、変数ｘｉが転送される。各計算ノードの共有メモリには、例えば、第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）、第２ベクトル（ｙ１，ｙ２，・・・，ｙＮ）のＬ個の変数、およびテンソルＪ（ｎ）の一部が保存される。そして、計算ノード間を接続する高速リンクでは、例えば、第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）が転送される。Ａｌｌｇａｔｈｅｒ関数を使う場合、各プロセッサで変数ｙｉを更新するために、第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）の全要素が必要となる。
　レプリカ交換法を使う場合、それぞれのレプリカにおける第１ベクトルおよび第２ベクトルの計算のために１または複数の計算ノードを割り当てることができる。この場合、レプリカ交換処理の実行タイミングでは、高速リンクを介して第１ベクトルおよび第２ベクトルのデータが計算ノード間で転送される。例えば、管理サーバ１の制御部１３が計算ノードの制御を行うことができる。なお、図１５に示したデータの配置および転送は一例にしかすぎない。ＰＣクラスタにおけるデータの配置方法、転送方法および並列化の実現方法については、特に問わない。
　また、ＧＰＵ（Ｇｒａｐｈｉｃｓ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）を使ってシミュレーテッド分岐アルゴリズムの計算を行ってもよい。
　図１６は、ＧＰＵを使った構成の例を概略的に示している。図１６には、互いに高速リンクで接続された複数のＧＰＵが示されている。それぞれのＧＰＵには、共有メモリにアクセス可能な複数のコアが搭載されている。また、図１６の構成例では、複数のＧＰＵが高速リンクを介して接続されており、ＧＰＵクラスタを形成している。例えば、ＧＰＵが図１のそれぞれの計算サーバに搭載されている場合、高速リンクは、ケーブル４ａ～４ｃおよびスイッチ５によって形成された計算サーバ間のインターコネクトに相当する。なお、図１６の構成例では、複数のＧＰＵが使われているが、ひとつのＧＰＵを使った場合にも、並列的な計算を実行することが可能である。すなわち、図１６のそれぞれのＧＰＵは、図１５のそれぞれの計算ノードに相当する計算を実行できる。すなわち、情報処理装置（計算サーバ）のプロセッサは、Ｇｒａｐｈｉｃｓ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ（ＧＰＵ）のコアであってもよい。
　ＧＰＵにおいて、変数ｘｉおよびｙｉ、ならびにテンソルＪ（ｎ）はデバイス変数として定義される。ＧＰＵは、変数ｙｉの更新に必要なテンソルＪ（ｎ）と第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）の積を、行列ベクトル積関数によって並列的に計算することができる。なお、行列とベクトルの積演算を繰り返し実行することにより、テンソルとベクトルの積を求めることができる。また、第１ベクトル（ｘ１，ｘ２，・・・，ｘＮ）の計算と、第２ベクトル（ｙ１，ｙ２，・・・，ｙＮ）のうち、積和演算以外の部分については、それぞれのスレッドにｉ番目の要素（ｘｉ，ｙｉ）の更新処理を実行させ、処理の並列化を実現することができる。
［組合せ最適化問題を解くための全体的な処理］
　以下では、シミュレーテッド分岐アルゴリズムを用いて組合せ最適化問題を解くために実行される全体的な処理を説明する。
　図１７のフローチャートは、組合せ最適化問題を解くために実行される全体的な処理の例を示している。以下では、図１７を参照しながら、処理を説明する。
　はじめに、組合せ最適化問題を定式化する（ステップＳ２０１）。そして、定式化された組合せ最適化問題をイジング問題（イジングモデルの形式）に変換する（ステップＳ２０２）。次に、イジングマシン（情報処理装置）によってイジング問題の解を計算する（ステップＳ２０３）。そして、計算された解を検証する（ステップＳ２０４）。例えば、ステップＳ２０４では、制約条件が満たされているか否かの確認が行われる。また、ステップＳ２０４でエネルギー関数（ハミルトニアン）の値を参照し、得られた解が最適解またはそれに近い近似解であるか否かの確認を行ってもよい。
　そして、ステップＳ２０４における検証結果または計算回数の少なくともいずれかに応じて再計算をするか否かを判定する（ステップＳ２０５）。再計算をすると判定された場合（ステップＳ２０５のＹＥＳ）、ステップＳ２０３およびＳ２０４の処理が再び実行される。一方、再計算をしないと判定された場合（ステップＳ２０５のＮＯ）、解の選択を行う（ステップＳ２０６）。例えば、ステップＳ２０６では、制約条件の充足またはエネルギー関数の値の少なくともいずれかに基づき選択を行うことができる。なお、複数の解が計算されていない場合には、ステップＳ２０６の処理をスキップしてもよい。最後に、選択した解を組合せ最適化問題の解に変換し、組合せ最適化問題の解を出力する（ステップＳ２０７）。
　上述で説明した情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラムを使うことにより、組合せ最適化問題の解を実用的な時間内で計算することが可能となる。これにより、組合せ最適化問題の求解がより容易となり、社会のイノベーションおよび科学技術の進歩を促進することが可能となる。
　なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組合せにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組合せてもよい。
１　管理サーバ
２　ネットワーク
３ａ、３ｂ、３ｃ　計算サーバ
４ａ、４ｂ、４ｃ　ケーブル
５　スイッチ
６　クライアント端末
１０　プロセッサ
１１　管理部
１２　変換部
１３　制御部
１４　記憶部
１４Ａ　問題データ
１４Ｂ　計算データ
１４Ｃ　管理プログラム
１４Ｄ　変換プログラム
１４Ｅ、３４Ｃ　制御プログラム
１５、３１　通信回路
１６　入力回路
１７　出力回路
１８　操作装置
１９　表示装置
２０　バス
３２　共有メモリ
３３Ａ、３３Ｂ、３３Ｃ、３３Ｄ　プロセッサ
３４　ストレージ
３４Ａ　計算データ
３４Ｂ　計算プログラム
３５　ホストバスアダプタ

Claims (15)

  1. 　それぞれが第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている、複数の演算回路と、
    　データ交換回路とを備え、
    　それぞれの前記演算回路は、前記第１変数を対応する第２変数に基づいて更新し、前記第１変数を第１係数で重み付けし対応する前記第２変数に加算し、複数の前記第１変数を用いて問題項を計算し、前記問題項を前記第２変数に加算するように構成されており、
    　それぞれの前記演算回路において前記第１係数に異なる値が設定されており、
    　前記データ交換回路は、前記演算回路間で前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されており、
    　情報処理装置。
  2. 　前記データ交換回路は、前記演算回路のペアごとに、条件判定を実行し、前記条件判定が肯定的であるときに前記演算回路間で計算に使われる前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されている、
    　請求項１に記載の情報処理装置。
  3. 　前記データ交換回路は、前記ペアに含まれる前記演算回路における前記第１ベクトル、前記第２ベクトルおよび前記第１係数に基づいて前記条件判定を実行するように構成されている、
    　請求項２に記載の情報処理装置。
  4. 　それぞれの前記演算回路で使われる前記第１係数は、互いに等比数列の関係となっている、
    　請求項１ないし３のいずれか一項に記載の情報処理装置。
  5. 　それぞれの前記演算回路は、前記第１ベクトルおよび前記第２ベクトルの更新回数に応じて、前記第１係数の値を周期的に変化させる、
    　請求項１ないし３のいずれか一項に記載の情報処理装置。
  6. 　前記データ交換回路は、タイミングによって異なる前記演算回路のペアを選択し、前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されている、
    　請求項１に記載の情報処理装置。
  7. 　前記データ交換回路は、複数の前記演算回路において、第１組合せの前記ペアと、第２組合せの前記ペアとを交互に選択し、前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されている、
    　請求項６に記載の情報処理装置。
  8. 　前記演算回路が計算する前記問題項は、イジングモデルに基づいている、
    　請求項１ないし７のいずれか一項に記載の情報処理装置。
  9. 　前記演算回路が計算する前記問題項は、多体相互作用を含んでいる、
    　請求項８に記載の情報処理装置。
  10. 　それぞれが第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている、複数の第１計算機と、
    　前記第１計算機間におけるデータの交換を実行するように構成された第２計算機とを備え、
    　それぞれの前記第１計算機は、第１変数を対応する第２変数に基づいて更新し、前記第１変数を第１係数で重み付けし対応する前記第２変数に加算し、複数の前記第１変数を用いて問題項を計算し、前記問題項を前記第２変数に加算するように構成されており、
    　前記第２計算機は、前記第１計算機間で前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されており、
    　それぞれの前記第１計算機において前記第１係数に異なる値が設定されている、
    　情報処理システム。
  11. 　前記第２計算機は、前記第１計算機のペアごとに、条件判定を実行し、前記条件判定が肯定的であるときに前記第１計算機間で計算に使われる前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するように構成されている、
    　請求項１０に記載の情報処理システム。
  12. 　それぞれが第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている、複数の第１計算機と、前記第１計算機間におけるデータの交換を実行するように構成された第２計算機とを使った情報処理方法であって、
    　それぞれの前記第１計算機が第１変数を対応する第２変数に基づいて更新するステップと、
    　それぞれの前記第１計算機が前記第１変数を第１係数で重み付けし対応する前記第２変数に加算するステップと、
    　それぞれの前記第１計算機が複数の前記第１変数を用いて問題項を計算するステップと、
    　それぞれの前記第１計算機が前記問題項を前記第２変数に加算するステップと、
    　前記第２計算機が前記第１計算機間で前記第１ベクトルおよび前記第２ベクトルの要素の交換、または前記第１係数の交換の少なくともいずれかを実行するステップとを含む、
    　情報処理方法。
  13. 　それぞれが第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新するように構成されている、複数の第１プロセッサと、前記第１プロセッサ間におけるデータの交換を実行するように構成された第２プロセッサとを使った情報処理方法であって、
    　それぞれの前記第１プロセッサが第１変数を対応する第２変数に基づいて更新するステップと、
    　それぞれの前記第１プロセッサが前記第１変数を第１係数で重み付けし対応する前記第２変数に加算するステップと、
    　それぞれの前記第１プロセッサが複数の前記第１変数を用いて問題項を計算するステップと、
    　それぞれの前記第１プロセッサが前記問題項を前記第２変数に加算するステップと、
    　前記第２プロセッサが前記第１プロセッサ間で前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するステップとを含む、
    　情報処理方法。
  14. 　複数の第１プロセッサのそれぞれに第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新させ、第２プロセッサが前記第１プロセッサ間におけるデータの交換を実行させるプログラムであって、
    　複数の第１プロセッサのそれぞれが第１変数を対応する第２変数に基づいて更新するステップと、
    　それぞれの前記第１プロセッサが前記第１変数を第１係数で重み付けし対応する前記第２変数に加算するステップと、
    　それぞれの前記第１プロセッサが複数の前記第１変数を用いて問題項を計算するステップと、
    　それぞれの前記第１プロセッサが前記問題項を前記第２変数に加算するステップと、
    　第２プロセッサが前記第１プロセッサ間で前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するステップとをコンピュータに実行させるプログラムを格納している、
    　非一時的なコンピュータ可読な記憶媒体。
  15. 　複数の第１プロセッサのそれぞれに第１変数を要素とする第１ベクトルおよび前記第１変数に対応する第２変数を要素とする第２ベクトルを繰り返し更新させ、第２プロセッサが前記第１プロセッサ間におけるデータの交換を実行させるプログラムであって、
    　複数の第１プロセッサのそれぞれが第変数を対応する第２変数に基づいて更新するステップと、
    　それぞれの前記第１プロセッサが前記第１変数を第１係数で重み付けし対応する前記第２変数に加算するステップと、
    　それぞれの前記第１プロセッサが複数の前記第１変数を用いて問題項を計算するステップと、
    　それぞれの前記第１プロセッサが前記問題項を前記第２変数に加算するステップと、
    　第２プロセッサが前記第１プロセッサ間で前記第１ベクトルおよび前記第２ベクトルの交換、または前記第１係数の交換の少なくともいずれかを実行するステップとを含む、
    　プログラム。
PCT/JP2020/014190 2019-03-28 2020-03-27 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム WO2020196872A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CA3135142A CA3135142A1 (en) 2019-03-28 2020-03-27 Information processing device, information processing system, information processing method, storage medium, and program
CN202080025407.3A CN113646787A (zh) 2019-03-28 2020-03-27 信息处理装置、信息处理系统、信息处理方法、存储介质及程序
JP2021509666A JP7421545B2 (ja) 2019-03-28 2020-03-27 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム
US17/487,179 US20220012017A1 (en) 2019-03-28 2021-09-28 Information processing device, information processing system, information processing method, and storage medium

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2019064699 2019-03-28
JP2019-064699 2019-03-28

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US17/487,179 Continuation US20220012017A1 (en) 2019-03-28 2021-09-28 Information processing device, information processing system, information processing method, and storage medium

Publications (1)

Publication Number Publication Date
WO2020196872A1 true WO2020196872A1 (ja) 2020-10-01

Family

ID=72608459

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/014190 WO2020196872A1 (ja) 2019-03-28 2020-03-27 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム

Country Status (5)

Country Link
US (1) US20220012017A1 (ja)
JP (1) JP7421545B2 (ja)
CN (1) CN113646787A (ja)
CA (1) CA3135142A1 (ja)
WO (1) WO2020196872A1 (ja)

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
HAYATO GOTO ET AL.: "Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems", SCIENCE ADVANCE, vol. 5, no. 4, 19 April 2019 (2019-04-19), pages 1 - 8, XP055715980, Retrieved from the Internet <URL:https://advances.sciencemag.org/content/5/4/eaav2372> *
HAYATO GOTO: "Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network", SCIENTIFIC REPORT, 22 February 2016 (2016-02-22), pages 1 - 8, XP055744351, Retrieved from the Internet <URL:https://www.nature.com/articles/srep21686> *
KOSUKE TATSUMURA ET AL.: "FPGA-Based Simulated Bifurcation Machine", 29TH INTERNATIONAL CONFERENCE ON FIELD PROGRAMABLE LOGIC AND APPLICATIONS (FPL), 2019, pages 59 - 66, XP033649156, Retrieved from the Internet <URL:https://ieeexplore.ieee.org/document/8892209> DOI: 10.1109/FPL.2019.00019 *

Also Published As

Publication number Publication date
US20220012017A1 (en) 2022-01-13
JPWO2020196872A1 (ja) 2020-10-01
CN113646787A (zh) 2021-11-12
CA3135142A1 (en) 2020-10-01
JP7421545B2 (ja) 2024-01-24

Similar Documents

Publication Publication Date Title
US11410070B2 (en) Syndrome data compression for quantum computing devices
WO2020196862A1 (ja) 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム
JP7421291B2 (ja) 情報処理装置、プログラム、情報処理方法、および電子回路
WO2019216277A1 (ja) 情報処理装置、演算装置、及び情報処理方法
JP7007520B6 (ja) 情報処理装置、演算装置、及び情報処理方法
Fang et al. RedSync: reducing synchronization bandwidth for distributed deep learning training system
Gyoten et al. Area efficient annealing processor for ising model without random number generator
US20220012307A1 (en) Information processing device, information processing system, information processing method, and storage medium
WO2020196883A1 (ja) 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム
JP7341965B2 (ja) 計算装置、計算方法およびプログラム
WO2020196872A1 (ja) 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム
Cardamone et al. Field‐programmable gate arrays and quantum Monte Carlo: Power efficient coprocessing for scalable high‐performance computing
WO2020196915A1 (ja) 情報処理装置、情報処理システム、情報処理方法、記憶媒体およびプログラム
JP7472062B2 (ja) 計算装置、計算方法およびプログラム
WO2022249785A1 (ja) 求解装置、求解方法およびプログラム
CN117313884A (zh) 量子电路处理方法、装置及电子设备
CN117313879A (zh) 量子电路处理方法、装置及电子设备

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20777978

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2021509666

Country of ref document: JP

Kind code of ref document: A

ENP Entry into the national phase

Ref document number: 3135142

Country of ref document: CA

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20777978

Country of ref document: EP

Kind code of ref document: A1