Quantum monte carlo calculation method, quantum monte carlo calculation device, computer program, and method for solving route optimization problem

Abstract

A quantum Monte Carlo calculation method according to the present invention comprises: a step for preparing a prescribed number of spin arrays for the Ising model; a step for setting an initial state for a minimum energy state pertaining to the Ising model, and an initial value for minimum energy; and a Monte Carlo calculation step for repeating a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until an end condition is satisfied. The Monte Carlo calculation step includes a step for generating, for each of the prescribed number of spin arrays, a new state candidate that is the spin state of the spin array to be processed with a change having been made in said spin state, and a step for performing a process of updating the state of the spin array to be processed and the minimum energy state on the basis of energy pertaining to the state candidate, energy pertaining to the spin array to be processed, and the minimum energy.

Description

Quantum Monte Carlo calculation method, quantum Monte Carlo calculation device, computer program, and method for solving path optimization problems

[0001] This disclosure relates to a quantum Monte Carlo calculation method, a quantum Monte Carlo calculation device, a computer program, and a method for solving a path optimization problem. This application claims priority to Japanese Application No. 2023-117400, filed on July 19, 2023, and incorporates the entire disclosure of said Japanese application by reference.

[0002] One type of quantum computing technology is called quantum annealing. Quantum annealing is specialized for solving combinatorial optimization problems. Some quantum annealing-based quantum computers have already been put into practical use. Quantum annealing is a metaheuristic method, and output of the optimal solution is not guaranteed. However, compared to quantum gate methods, which are other quantum computing methods, at the current technological level, it is expected to be able to output the optimal solution with a higher probability.

[0003] Quantum annealing uses a model known as the Ising model. This Ising model is described as an undirected graph containing multiple nodes, each with a spin, and edges between these nodes. Quantum annealing determines the spin configuration that gives the Ising model the minimum energy, assuming that the spins of these nodes are quantum mechanically superposed. The path optimization problem is a type of combinatorial optimization problem that is difficult to solve, but it is known that a candidate solution can be obtained by determining the spin configuration in the minimum-energy state of such an Ising model.

[0004] There is a method called the quantum Monte Carlo method for simulating this quantum annealing using a classical computer. For the Ising model in quantum annealing, a method called Suzuki-Trotter decomposition can be used. By this method, a concept called imaginary time is introduced, and for example, an Ising model represented by a two-dimensional graph can be analyzed by a three-dimensional classical model with one higher dimension. Each of these layers is called a Trotter layer. By applying the Monte Carlo method to these Trotter layers, an estimation is made of the spin configuration that gives the minimum energy of the original Ising model. Such a technique is disclosed in Non-Patent Document 1 shown below.

[0005] Zhang, Hengshuo, et al., "Experimental Evaluation of the Path-Integral Monte Carlo Method for Restricted Boltzmann Machines", The 104th Symposium on Fundamental Problems of Artificial Intelligence, The Japanese Society for Artificial Intelligence, Session: 104, Date of Meeting: August 8 - 9, 2017, [online], [searched on June 1, 2023], Internet <URL: https: / / www.jstage.jst.go.jp / article / jsaifpai / 104 / 0 / 104_02 / _article / -char / ja / >Masuda, Kenichi, et al., "Optimization of Distribution Planning by Quantum Computing", Sumitomo Electric Technical Review, Vol. 202, January 2023, p. 1 - 4Hasitha Muthumala Widyasooriyana, et al., "A GPU-based quantum annealing simulator for fully-connected ising models utilizing spatial and temporal parallelism", [online], IEEE Access, Vol. 8, 9057502, 2020, p. 67929 - 67939. Internet <url: https: / / doi.org / 10.1109 / ACCESS.2020.<2985699>

[0006] A quantum Monte Carlo calculation method according to a first aspect of this disclosure is a quantum Monte Carlo calculation method for calculating a minimum energy state that provides minimum energy, which is the minimum value of energy in an Ising model, and includes: a step in which a computer prepares a predetermined number of spin arrangements for the Ising model; a step in which a computer sets an initial state of the minimum energy state for the Ising model and an initial value of the minimum energy; and a Monte Carlo calculation step in which the computer repeats a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until a termination condition is met, wherein the Monte Carlo calculation step includes: a step in which, for each of the predetermined number of spin arrangements, the computer generates a new state candidate in which the spin state of the spin arrangement to be processed has changed; and a step in which the computer performs processing to update the state of the spin arrangement to be processed and the minimum energy state based on the energy for the state candidate, the energy for the spin arrangement to be processed, and the minimum energy obtained up to now.

[0007] FIG. 1 is a schematic diagram of a two-dimensional Ising model. FIG. 2 is a conceptual diagram showing the state in the imaginary time direction after Suzuki-Trotter decomposition of a two-dimensional Ising model. FIG. 3 is a flowchart showing the structure of a main routine of a program implementing a route optimization problem solving device employing a quantum Monte Carlo calculation method according to an embodiment of the present disclosure. FIG. 4 is a flowchart showing the control structure of a program implementing the quantum Monte Carlo calculation of the program shown in FIG. 3. FIG. 5 is a graph showing the control of temperature and transverse magnetic field control in the quantum Monte Carlo calculation method according to the first embodiment of the present disclosure. FIG. 6 is a flowchart showing the control structure of a program implementing the steps of updating the minimum energy and spin state of the Ising model, among the programs shown in FIG. 4. FIG. 7 is a diagram showing, in tabular form, the results of solving a benchmark problem of a max-cut problem described in Non-Patent Document 2 using the quantum Monte Carlo method (quantum annealing simulator). FIG. 8 is a diagram showing, in tabular form, the results of solving the same benchmark problem shown in FIG. 7 according to the first embodiment of the present disclosure. FIG. 9 is a graph showing the convergence of a max-cut problem for benchmark problem dataset G34 according to the first embodiment of the present disclosure, compared with the results obtained using a conventional simulated annealing method. FIG. 10 is a block diagram showing the hardware configuration of a computer system for realizing one embodiment of this disclosure. FIG. 11 is a graph showing the control of temperature and transverse magnetic field in a second embodiment of this disclosure. FIG. 12 is a flowchart showing the control structure of a program for realizing quantum Monte Carlo calculation according to the second embodiment of this disclosure. FIG. 13 is a flowchart showing the control structure of a program for realizing quantum energy and state update processing, among the programs whose control structure is shown in FIG. 12. FIG. 14 is a graph showing the convergence of a max-cut problem for dataset G9 according to the second embodiment of this disclosure, compared with the results according to the first embodiment. FIG. 15 is a graph showing the convergence of a max-cut problem for dataset G20 according to the second embodiment of this disclosure, compared with the results according to the first embodiment.16 is a graph showing the convergence of the max-cut problem for data set G34 according to the second embodiment of this disclosure, compared with the results obtained by a conventional simulated annealing (SA) method. FIG. 17 is a graph showing the convergence of the max-cut problem for data set G34 according to the second embodiment of this disclosure, compared with the results obtained by the first embodiment of this disclosure.

[0008] [Problems to be Solved by This Disclosure] However, conventional quantum Monte Carlo calculation methods have had problems such as requiring a long time for calculation and poor convergence.

[0009] Therefore, an object of this disclosure is to provide a quantum Monte Carlo calculation method, a quantum Monte Carlo calculation device, a computer program, and a method for solving a path optimization problem that can shorten the calculation time.

[0010] Effect of the Disclosure As described above, according to the disclosure, it is possible to provide a quantum Monte Carlo calculation method, a quantum Monte Carlo calculation device, a computer program, and a method for solving a route optimization problem that can reduce calculation time.

[0011] [Description of the embodiments of the present disclosure] In the following description and drawings, the same parts are denoted by the same reference numerals. Therefore, detailed description thereof will not be repeated. Note that at least some of the embodiments described below may be combined in any manner.

[0012] (1) A quantum Monte Carlo calculation method according to a first aspect of this disclosure is a quantum Monte Carlo calculation method for calculating a minimum energy state that provides minimum energy, which is the minimum value of energy in an Ising model, the method including: a step in which a computer prepares a predetermined number of spin arrangements for the Ising model; a step in which a computer sets an initial state of the minimum energy state for the Ising model and an initial value of the minimum energy; and a Monte Carlo calculation step in which the computer repeats a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until a termination condition is met, wherein the Monte Carlo calculation step includes a step in which, for each of the predetermined number of spin arrangements, the computer generates a new state candidate in which the spin state of the spin arrangement to be processed has changed; and a step in which the computer performs processing to update the state of the spin arrangement to be processed and the minimum energy state based on the energy for the state candidate, the energy for the spin arrangement to be processed, and the minimum energy obtained up to now.

[0013] With this configuration, the Trotter layer state and the minimum energy state are updated whenever the energy for a certain spin configuration and the minimum energy satisfy the conditions. Compared to calculating the energies for all Trotter layers and then attempting to update the minimum energy state, this reduces the amount of calculation and enables the minimum energy state to be calculated quickly and with good convergence.

[0014] (2) In the above (1), the energy of the Ising model is a sum of a first energy due to interactions within the spin array to be processed and a second energy due to interactions between the spin array to be processed and the spin array that gives the minimum energy, and the step of performing the update process includes a first determination step in which a computer determines whether or not a difference between the first energy related to the state candidate and the first energy related to the spin array to be processed is smaller than 0, and a second determination step in which a computer updates the spin array to be processed based on the first energy related to the state candidate and the minimum energy obtained up to the present time, in response to a positive determination in the first determination step. The method may include a first selective processing step in which a process of updating the state of the array using the state candidate and updating the minimum energy state and the minimum energy using the state candidate and the first energy related to the state candidate, respectively, and a process of updating the state of the spin array to be processed using the state candidate, and a second selective processing step in which, in response to a negative judgment in the first judgment step, the computer selectively executes a process of updating or not updating the state of the spin array to be processed using the state candidate, depending on whether the difference between the energy related to the state candidate and the energy related to the spin array to be processed is less than 0.

[0015] With this configuration, the update processes for each of the predetermined number of spin configurations can be performed independently of each other, and the minimum energy state obtained in each update process can be used in other update processes. Therefore, the minimum energy state can be updated more efficiently than when the minimum energy update is performed after the energy calculations for all spin configurations are completed. As a result, the calculation of the minimum energy state can be performed quickly, in a short time, and with good convergence.

[0016] (3) In the above (2), a contribution of the first energy to the sum of the first energy and the second energy may be greater than a contribution of the second energy.

[0017] By giving more weight to the first energy than the second energy while appropriately selecting the contribution of the second energy, higher convergence is achieved.

[0018] (4) In the above (3), the contribution of the first energy may be set to 1, and the contribution of the second energy may be set to 0.

[0019] By considering only the first energy, the convergence is highest and the minimum energy state can be calculated quickly.

[0020] (5) In any one of (2) to (4) above, the first selective processing step may include: a second determination step in which, in response to a positive determination in the first determination step, the computer determines whether the first energy associated with the state candidate is smaller than the first energy associated with the spin array to be processed; a first update step in which, in response to a positive determination in the second determination step, the computer updates the state of the spin array to be processed with the state candidate and updates the minimum energy state and the minimum energy with the state candidate and the first energy associated with the state candidate, respectively; and a second update step in which, in response to a negative determination in the second determination step, the computer updates the state of the spin array to be processed with the state candidate.

[0021] With this configuration, the process of updating both the minimum energy state and the spin arrangement state based on the change in the first energy due to the interaction between the nodes of each Trotter layer, and the process of updating the states of each Trotter layer are appropriately selected and executed so that the energy decreases toward the minimum energy state, thereby enabling the calculation of the minimum energy state to be performed quickly, in a short time, and with good convergence.

[0022] (6) In (5) above, the second selective processing step may include a third judgment step in which, in response to a negative judgment in the first judgment step, the computer judges whether the difference between the energy related to the state candidate and the energy related to the spin array to be processed is greater than or equal to 0, and a third update step in which, in response to a positive judgment in the third judgment step, the computer performs processing to update the state of the spin array to be processed using the state candidate.

[0023] This configuration allows the spin ordering of the Trotter layers to be updated to align the spin orderings not only when the energy is decreasing but also when the energy is increasing, thereby eliminating the possibility that the energy of the entire system will become stuck at a local minimum.

[0024] (7) In (6) above, the third update step may include a step in which a computer updates the state of the spin array to be processed with the state candidate according to a probability determined by the difference between the first energy related to the state candidate and the first energy related to the spin array to be processed, i.e., a step in which the computer updates the state of the spin array using the Metropolis method.

[0025] With this configuration, even if the first energy increases in the energy calculation of each Trotter layer, the spin configuration state of the Trotter layer is updated with a candidate state with a certain probability as long as the difference in the total energy is equal to or less than zero. As a result, the possibility that the energy of the entire system will remain at a local minimum can be reduced while maintaining good convergence. Furthermore, the probability is gradually reduced at this time. As a result, updates to spin configurations with increasing energy can be reduced as the calculation proceeds, allowing the solution to converge.

[0026] (8) In any one of the above (1) to (4), the Monte Carlo calculation step is performed in parallel for the predetermined number of spin arrangements.

[0027] The first energy is calculated only for the spin configuration of each Trotter layer that is the target of processing, while the second energy is calculated for both the spin configuration of the Trotter layer that is the target of processing and the spin configuration of the candidate state that provides the minimum energy value. Therefore, when updating the spin configuration of each Trotter layer, it is not necessary to consider the spin configurations of other Trotter layers. Therefore, the Monte Carlo calculation steps can be performed in parallel. As a result, not only is the sum calculation required in the prior art unnecessary, but by parallelizing the update of the spin configurations of each Trotter layer, the Monte Carlo calculation steps can be performed in parallel, allowing for faster solution retrieval.

[0028] (9) In any one of (1) to (4) above, the step of generating the new state candidate may include a step in which a computer randomly selects one of the nodes of the spin array to be processed and inverts the spin of that node.

[0029] This configuration can reduce the possibility that the energy of the entire system will stagnate at a local minimum value while maintaining good convergence.

[0030] (10) In the above (1), the energy of the Ising model may be the sum of a first energy due to interactions within the spin array to be processed and a second energy due to interactions between the spin array to be processed and the spin array that gives the minimum energy obtained up to the present time, and the step of performing the update processing may include an update step in which a computer executes a process of updating the minimum energy with the energy of the spin array to be processed in response to the first energy for the state candidate being smaller than the minimum energy obtained up to the present time, a process of updating the state of the spin array to be processed with the state candidate, and a process of updating the spin array with the state candidate in accordance with a probability determined by a difference between the energy for the state candidate and the energy for the spin array to be processed.

[0031] With this configuration, the update processes for each of the predetermined number of spin configurations can be performed independently of each other, and the minimum energy state obtained in each update process can be used in other update processes. This allows the minimum energy state to be updated in a shorter time than if the minimum energy update was performed after the energy calculations for all spin configurations were completed. As a result, the calculation of the minimum energy state can be performed quickly and with good convergence.

[0032] (11) A quantum Monte Carlo calculation device according to a second aspect of this disclosure is a quantum Monte Carlo calculation device for calculating a minimum energy state that gives a minimum energy in an Ising model, the quantum Monte Carlo calculation device including: a Trotter layer preparation unit that prepares a predetermined number of spin arrangements for the Ising model; an initial state setting unit that sets an initial state of the minimum energy state for the Ising model and an initial value of the minimum energy; and a Monte Carlo calculation unit that repeats a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until a termination condition is met, wherein the Monte Carlo calculation unit includes: a state candidate generation unit that generates, for each of the predetermined number of spin arrangements, a new state candidate in which the spin state of the spin arrangement to be processed has changed; and an update processing unit that performs, for each of the predetermined number of spin arrangements, an update process of the state of the spin arrangement to be processed and the minimum energy state based on the energy related to the state candidate, the energy related to the spin arrangement to be processed, and the minimum energy obtained up to now.

[0033] With this configuration, if the energy of a certain spin configuration and the minimum energy satisfy the conditions, the Trotter layer state and the minimum energy state are updated. Compared to calculating the sum of the energies of all Trotter layers and then attempting to update the minimum energy state, this reduces the amount of calculation, allowing the minimum energy state to be calculated quickly and with good convergence.

[0034] (12) A computer program according to a third aspect of this disclosure is a computer program that causes a computer to execute a quantum Monte Carlo calculation method for calculating a minimum energy state that gives a minimum energy in an Ising model, the quantum Monte Carlo calculation method including: a step in which a computer prepares a predetermined number of spin arrangements for the Ising model; a step in which a computer sets an initial state of the minimum energy state for the Ising model and an initial value of the minimum energy; and a Monte Carlo calculation step in which a computer repeats a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until a termination condition is met, the Monte Carlo calculation step including: a step in which, for each of the predetermined number of spin arrangements, a computer generates a new state candidate in which the spin state of the spin arrangement to be processed has changed; and a step in which a computer performs processing to update the state of the spin arrangement to be processed and the minimum energy state based on the energy for the state candidate, the energy for the spin arrangement to be processed, and the minimum energy obtained up to now.

[0035] With this configuration, if the energy of a certain spin configuration and the minimum energy satisfy the conditions, the Trotter layer state and the minimum energy state are updated. Compared to calculating the sum of the energies of all Trotter layers and then attempting to update the minimum energy state, this reduces the amount of calculation, allowing the minimum energy state to be calculated quickly and with good convergence.

[0036] (13) A method for solving a route optimization problem according to a fourth aspect of this disclosure includes the steps of: converting, by a computer, a route optimization problem formulated in a predetermined format into a max-cut problem and an Ising model equivalent to the max-cut problem; determining, by the computer, a solution to the Ising model by applying the quantum Monte Carlo calculation method described in (1) above to the Ising model; and converting, by the computer, the solution into a solution to the route optimization problem.

[0037] This configuration allows the route optimization problem to be solved in a shorter time than conventional methods, with high convergence. As a result, it is more likely that packages and other items can be delivered in a timely and efficient manner along the optimal route, whereas conventional methods have not been able to deliver packages and other items in a timely manner along the optimal route.

[0038] [Details of the embodiments of the present disclosure] Specific examples of a quantum Monte Carlo calculation method, a quantum Monte Carlo calculation device, a computer program, and a method for solving a path optimization problem according to embodiments of the present disclosure will be described below with reference to the drawings. Note that the present disclosure is not limited to these examples, but is defined by the scope of the claims, and is intended to include all modifications within the meaning and scope equivalent to the claims.

[0039] 1. First Embodiment 1.1 Ising Model Referring to FIG. 1, a two-dimensional Ising model 50 includes a plurality of nodes 60, 62, 64, etc., as described above, and a plurality of edges connecting these nodes. Spins 70, 72, etc. are assigned to each of the nodes 60, 62, 64, etc. In classical physics, spins have values ​​of (1, -1), which respectively represent an up-point and a down-point. On the other hand, in quantum mechanics, spins are treated as a superposition of the two.

[0040] FIG. 2 shows the state of this Ising model 50 after Suzuki-Trotter decomposition. Referring to FIG. 2, the Ising model 90 after Suzuki-Trotter decomposition consists of multiple Trotter layers 100, 102, ..., 104. These Trotter layers 100, 102, ..., 104 are spin arrangements decomposed in the direction of an axis called the Trotter axis (also called the imaginary time axis). In conventional quantum Monte Carlo methods, the minimum energy state is found by calculating the sum of the energies of these Trotter layers 100, 102, ..., 104 using a classical computer.

[0041] First, the Hamiltonian H representing the energy of the system of the Ising model 50 is expressed by the following equation.

[0042] In this formula, J ｉｊindicates the magnitude of the interaction between adjacent spins in the spin arrangement set as the problem. ｉ ｘ and σ ｉ ｚ are the Pauli spin matrices, and are expressed by the following matrices, respectively:

[0043] Γ represents the magnitude of the transverse magnetic field applied to the Ising model.

[0044] The partition function Z can be calculated from the Hamiltonian H using the following equation (3): Ｂ T(k Ｂ is the Boltzmann constant, and T is temperature. In SA, β changes in proportion to the inverse of the temperature. Therefore, arbitrary control is allowed in the calculation. In the first embodiment, this β is changed. In the second embodiment described later, this β is kept constant.

[0045] The objective function is to find the energy E that maximizes this partition function Z, that is, to find the spin arrangement that minimizes the energy.

[0046] By applying Suzuki-Trotter decomposition to the above equation (1), it can be expressed by a three-dimensional classical expression in which an imaginary time axis is introduced.

[0047] In this formula, P represents the number of Trotter layers, and s ｋ ｉ represents the i-th spin (1 or −1) of the k-th Trotter layer. ⊥ represents the magnitude of the interaction between the spins of different (adjacent) Trotter layers, and T represents the temperature during annealing.

[0048] As a result, as described above, it comes down to finding the state that gives the minimum value of energy E in the Ising model, that is, the spin arrangement.

[0049] 1.2 Calculation Method In a typical quantum Monte Carlo calculation that faithfully calculates this equation, it is necessary to calculate the sum of the energies of each Trotter layer, as expressed by the following equation:

[0050] When calculating the above equation, the prior art required the following method to obtain the minimum energy value. Specifically, the spin configurations of all of the Trotter layers 100, 102, ..., 104 shown in FIG. 2 are updated independently. Furthermore, the results of adding these up are used to update the spin state and the minimum energy value of the Ising model 50. This calculation requires the calculation shown in equation (7) to be performed for all of the Trotter layers. Therefore, the prior art had problems in that each iteration of the calculation took a long time and convergence was poor.

[0051] In this first embodiment, this problem is solved as follows: argmin is calculated first, and the process of calculating the sum of the energies of all the Trotter layers is omitted.

[0052] In this formula, |α Ｋ 〉 is a new spin configuration candidate. This new spin configuration candidate is the Kth spin configuration |α Ｋ 〉', a predetermined number of spins (one in this embodiment) are randomly selected and inverted. ｍｉｎ 〉 is the spin configuration that gives the minimum energy at the current time among all P spin configurations of the Trotter layers, and σ ｚ ｍｉｎ，ｉ is the Pauli matrix acting on the spin configuration that gives the minimum energy. Therefore, E1 represents the energy given by the interaction between spins in the Kth Trotter layer. Also, E2 represents the energy given by the interaction between the spin configuration of the Kth Trotter layer and the spin configuration that currently gives the minimum energy.

[0053] According to this formula, the minimum value is not found by calculating the energy for all the Trotter layers. In this formula, the spin arrangement |α ｍｉｎ > is found, that value is used immediately in the calculation of other spin configurations. In this way, calculations for each Trotter layer can be performed in parallel, and when the minimum value is found, that value can be immediately reflected in the calculation for each Trotter layer. Therefore, the calculation method shown in this embodiment reduces the amount of calculation required to obtain the spin configuration that gives the minimum energy. As a result, a solution can be found quickly. In addition, the convergence to the spin configuration that gives the minimum energy is improved.

[0054] 1.3. Configuration This quantum Monte Carlo method can be applied to a variety of fields. One example is the optimization of parcel delivery plans in the logistics field. Delivery plan optimization is a type of combinatorial optimization problem, but as the number of vehicles used and the number of delivery destinations increase, the amount of calculations increases explosively, making it difficult to find an exact solution. The quantum Monte Carlo method described above has the potential to significantly reduce the time required to solve delivery plan optimization problems.

[0055] As described in Non-Patent Document 2, the above-described problem regarding delivery plan optimization can be formulated using a technique called QUBO (Quadratic Unconstrained Binary Optimization). It is also known that QUBO can be converted into a max-cut problem for a graph. The max-cut problem is a problem of finding a way to divide a graph into two subgraphs that maximizes (optimizes) the number of edges to be cut (when weights are assigned to the edges, the sum of the weights of the edges to be cut). The max-cut problem is known to be equivalent to the Ising model. Therefore, the above-described quantum Monte Carlo calculation method can be applied to the max-cut problem.

[0056] 3 is a flowchart showing the control structure of the main routine of a delivery plan optimization program to which the quantum Monte Carlo calculation method according to the first embodiment is applied. Referring to FIG. 3, this program includes step 150 of reading a QUBO formulation file that records a delivery plan optimization problem formulated in the QUBO format, and step 152 of converting the contents of the formulation file read in step 150 into a max-cut problem and an Ising model equivalent to this max-cut problem. This program further includes step 154 ​​of determining the spin configuration of the minimum energy state of the Ising model by applying the quantum Monte Carlo calculation method to the Ising model obtained in step 152, and step 156 of outputting the spin configuration calculated in step 154 ​​as a solution to the delivery plan optimization problem formulated as an initial QUBO problem and terminating the process.

[0057] Referring to FIG. 4, step 154 ​​shown in FIG. 3 performs Suzuki-Trotter decomposition on the Ising model corresponding to the max-cut problem to obtain initial spin arrangements (|α Ｋ 〉' (where K=1 to P)) are randomly generated and prepared in step 200, and the above |α Ｋ and step 204 for setting initial values ​​of the temperature (T) and transverse magnetic field strength (Γ) used as conditions for energy calculation in the following iterations.

[0058] Step 154 ​​further includes step 206, which repeats step 208, a process of calculating the minimum energy for all Trotter layers while changing the temperature T and transverse magnetic field Γ used as conditions for energy calculation, a predetermined number of times (i.e., until a termination condition is met), and step 210, which inversely transforms the spin arrangement that gives the minimum energy of the Ising model obtained at the end of step 208 into a solution to the delivery plan optimization problem formulated as QUBO, and outputs the solution, and terminates the process. In this first embodiment, as shown by the straight line 280 in Figure 5, the temperature T and transverse magnetic field Γ are each linearly updated, for example, each time step 208 is executed, so that they become very small values ​​from their predetermined initial values ​​after the above-mentioned predetermined number of repetitions. Furthermore, J obtained by equation (4) ⊥ The value of J is calculated by dividing the value by the temperature (T) at that time. ⊥ This kind of control is difficult for current quantum computers, but it is computationally possible. ⊥ 5 is obtained as a curve 282. As a result, as shown by the curve 282, the J ⊥ The contribution of

[0059] Step 208 includes step 230, which repeats step 232, which performs the minimum energy and spin alignment update, for the number of Trotter layers while changing the spin alignment of the Trotter layer of interest, and step 234, which updates the values ​​of temperature (T) and transverse magnetic field (Γ) after step 230, and terminates execution of step 208.

[0060] Step 232 includes step 250, in which a predetermined number (one in this embodiment) of spin arrangements in the target Trotter layer are randomly selected and inverted to generate new spin arrangement candidates (state candidates); step 252, in which energies E1 and E2 are calculated for the new state candidates; step 254, in which differences ΔE1 and ΔE2 are calculated between the new energies E1 and E2 calculated in step 252 and the energies E1 and E2 calculated from the current spin arrangements, respectively; and step 256, in which the minimum energy and the spin arrangement that gives the minimum energy are updated using the energies E1 and E2 calculated in step 252, ΔE1 and ΔE2 calculated in step 254, and the minimum energy calculated in the iterations up to the current time, thereby ending step 232.

[0061] At this time, if the minimum energy and the spin array that gives the minimum energy are updated in step 256, the spin array that gives the minimum energy may be used in step 256. There is no need to wait for the completion of step 232. Step 230 repeatedly executes step 232. However, in this embodiment, as described above, the control structure is such that the spin array that gives the minimum energy is updated without waiting for the completion of step 232. Therefore, if the iterative processing of step 232 is performed in parallel and a spin array that gives the minimum energy is obtained in any of the processes, that spin array can be used in all other processes to perform subsequent processing. In other words, this iterative processing can be performed in parallel on a GPU or the like, thereby shortening the time until a solution is obtained.

[0062] Steps 252, 254, and 256 update the Trotter layer state and the minimum energy state for each of a predetermined number of spin arrangements based on the energy associated with the candidate state, the energy associated with the spin arrangement being processed, and the minimum energy.

[0063] 6 shows the details of step 250 in Fig. 4. Referring to Fig. 6, step 250 includes step 300, which branches the flow of control according to whether ΔE1 is equal to or less than 0, and step 300, which branches the flow of control according to whether ΔE1 is equal to or less than 0, and step 300, which branches the flow of control according to whether ΔE1 is equal to or less than 0, in response to the determination in step 300 being affirmative. ｍｉｎ Step 250 further includes step 302 of branching the flow of control depending on whether the value is smaller than the minimum energy E. In this specification, "smaller or not" in relation to the magnitude of a numerical value may be interpreted as meaning either "equal to or less than the minimum energy E" or "less than the minimum energy E". In the case of FIG. 4, "smaller or not" means "equal to or less than the minimum energy E". Step 250 further includes step 302 of branching the flow of control depending on whether the value is smaller than the minimum energy E. ｍｉｎ is updated with the new energy E1, and the minimum energy E ｍｉｎ The spin configuration that gives the spin configuration candidate state |α Ｋ and in response to a negative determination in step 302, update the current spin arrangement |α Ｋ 〉' is a candidate state of spin arrangement |α Ｋ and step 306 of updating only to >. After steps 304 and 306, execution of step 250 ends.

[0064] Note that in step 300, a determination is made based only on ΔE1. ​​This is because, as the quantum Monte Carlo calculation progresses, E2 ultimately approaches zero, and the energy of the Ising model can be described by E1 alone. However, this disclosure is not limited to such an embodiment. For example, in step 300, the control flow may be branched depending on whether the value of the linear sum of ΔE1 and ΔE2 (a × ΔE1 + b × ΔE2) is equal to or less than 0. The above-described step 300 is an example in which a = 1 and b = 0.

[0065] In this way, the coefficients a and b in the above linear sum are determined to appropriately reflect the contribution of ΔE2 in the above calculation. In other words, these values ​​can also be determined so that a+b=1. In this case, the absolute values ​​of the coefficients a and b are not significant, but rather the ratio between the two is significant. In this way, the values ​​of (a, b) can be a=b=1, a+b>1, or a+b<1. Furthermore, a×ΔE １ ｎ + b × ΔE ｍ The form may be arbitrarily selected depending on the Ising model to be solved.

[0066] Step 250 further includes step 308, in response to the determination in step 300 being negative, of branching the flow of control according to whether or not ΔE1+ΔE2 (=difference of E1+difference of E2) is 0 or less; step 310, in response to the determination in step 308 being positive, of sampling a random number ξ from a uniform distribution of [0, 1]; step 312, in response to the determination in step 312 being positive, of branching the flow of control according to whether or not the value ξ sampled in step 310 is exp(−λΔE1 / T) or less; Ｋ 〉' as a state candidate |α Ｋ >, and terminating step 314. If the determination in step 308 is negative, or if the determination in step 312 is negative, the state change is not made and execution of step 250 terminates.

[0067] By doing this, even if the energy E1 increases due to the state candidate, the spin configuration of the Trotter layer is replaced with the state candidate with a certain probability as long as it approaches the spin configuration that minimizes the energy. As a result, the probability that the energy will remain at a local minimum can be reduced while maintaining convergence. Note that λ in the above equation can be adjusted appropriately according to parameters such as the number of Trotter layers.

[0068] As will be described later, the quantum Monte Carlo calculation method described above requires less time and obtains the best solution with a higher probability in benchmark max-cut problems compared to conventional methods. The reason for this is that, in the state update process for spin configurations, if a spin configuration that gives the minimum energy is found, that value is immediately used in subsequent calculations. This effectively reduces the number of iterative loops by one, compared to the conventional method of selecting the spin configuration with the minimum energy after completing calculations for all spin configurations, thereby reducing the amount of calculation.

[0069] The results of a benchmark experiment comparing the quantum Monte Carlo calculation method according to this embodiment with a conventional quantum Monte Carlo calculation method are shown below. For this experiment, a dataset called Gset, which is often used in this field as a benchmark for max-cut problems, was used.

[0070] FIG. 7 shows a table described in Non-Patent Document 3 as a result of applying the method described in Non-Patent Document 3 to datasets G9, G13, G18, G19, G20, and G21 of the benchmark datasets of Gset as a conventional technique. For example, dataset G9 is a max-cut problem for a graph containing 800 nodes. The number of sequence candidates that can be the optimal solution to this problem is 2 ８００ ≒10 ２４０ Therefore, it is extremely difficult to find an optimal solution for dataset G9. The same is true for the other datasets.

[0071] 7, the first column of the table shows the name of the benchmark dataset, and the second, third, and fourth columns show the number of nodes, the number of edges, and the numerical value obtained by the best known cut of the graph represented by the corresponding dataset, respectively.

[0072] The fifth column of this table shows the max cut values ​​obtained by 1000 MC (Monte Carlo) steps, as described in Non-Patent Document 3. The sixth and seventh columns show the best and average values ​​obtained by 1000 MC steps using a single-core CPU (Central Processing Unit). The eighth and ninth columns show the best and average values ​​obtained by 1000 MC steps when parallel calculations using a GPU (Graphics Processing Unit) are introduced, as described in Non-Patent Document 3.

[0073] FIG. 8 shows the results of applying the quantum Monte Carlo calculation method according to the embodiment to the benchmark data set using a quad-core CPU. The "temperature" shown in FIG. 8 is the initial value in the quantum Monte Carlo calculation method. The results shown in FIG. 7 were obtained by updating multiple spins at once and stopping the calculation when 1,000 MC steps were reached. On the other hand, the results shown in FIG. 8 were obtained by updating only one spin until the energy no longer decreased (approximately 500,000 MC steps), rather than stopping the calculation at a certain number of MC steps. In other words, the calculation conditions for the two methods are different, and the results shown in FIG. 8 are considered to be based on more favorable calculation conditions. Therefore, a simple comparison between the two is difficult.

[0074] However, there is still a significant difference between the results, as described below. For example, comparing Figures 7 and 8 for dataset G9 reveals the following facts. First, as shown in Figure 8, in this embodiment, the best known maxcut value (2054) was reached in 55 out of 100 trials. This result represents a significant improvement compared to the results shown in Figure 7. This improvement can also be seen from the number of datasets that reached the best known maxcut value. In Figure 7, only two datasets (G13 and G20) reached the best known value. In contrast, in the results shown in Figure 8, the best known value was reached for all datasets shown. Furthermore, when looking at the average values, Figure 8 exceeds Figure 7 except for G13. These results demonstrate that the quantum Monte Carlo calculation method according to the above embodiment can efficiently solve combinatorial optimization problems (e.g., delivery plan optimization problems) that were previously difficult to solve.

[0075] The results shown in FIG. 7 were obtained using the calculation method shown in Equation (7). On the other hand, the results shown in FIG. 8 were obtained using the calculation method shown in Equation (8). That is, the former calculates the sum of the energies and then finds the minimum energy state, while the latter calculates the minimum energy state and then finds the sum of the energies. There is no theoretical basis for whether correct results can be obtained by changing the calculation procedure in this way. However, the results shown in FIG. 8 show that sufficiently correct results can be obtained quickly by using the calculation method shown in Equation (8).

[0076] 9 shows the convergence of G34 according to the first embodiment (curve 330) and the convergence according to the conventional SA method (curve 332). As shown in FIG. 9, the first embodiment achieves faster convergence than the SA method. Moreover, as described above, the first embodiment achieves a higher rate of reaching the best known numerical value than the SA method.

[0077] 10 shows an example of the hardware configuration of a computer that executes a computer program for implementing the quantum Monte Carlo calculation method according to the above embodiment. The configuration shown in FIG. 10 is merely an example, and it goes without saying that the quantum Monte Carlo calculation method can also be implemented by a computer different from that shown in FIG. 10.

[0078] 10, this computer system 350 includes a computer 360, and a keyboard 366, a mouse 364, and a monitor 362 for interacting with a user, all of which are connected to the computer 360. Of course, these are just one example of a configuration for when user interaction is required, and any general hardware and software that can be used for user interaction (e.g., a touch panel, voice input, or a general pointing device) can be used.

[0079] The computer 360 includes a bus 410, a CPU 400, a GPU 402, a ROM (Read-Only Memory) 404 that stores a boot-up program for the computer 360, and the like, all of which are connected to the bus 410, a RAM (Random Access Memory) 406 that stores instructions constituting the program, system programs, working data, and the like when the program is executed, and an SSD (Solid State Drive) 408 that is non-volatile memory connected to the bus 410. The SSD 408 is used to store programs executed by the CPU 400 and GPU 402, as well as data used by the programs executed by the CPU 400 and GPU 402. These programs include a program for implementing the quantum Monte Carlo calculation method described above.

[0080] Computer 360 further includes a network I / F (Interface) 418 that provides connection to network 370, and a USB (Universal Serial Bus) port 416 to which a USB memory 368 can be attached and which provides communication between USB memory 368 and each component within computer 360. In this embodiment, network 370 is the Internet.

[0081] The computer 360 further includes an input / output I / F 414 connected to the bus 410 for inputting input from the keyboard 366 and the mouse 364 to the CPU 400, and a display control unit 412 connected to the bus 410 for controlling the display on the monitor 362 in accordance with control from the CPU 400.

[0082] A computer program for operating this computer system 350 to realize each function described below is stored in a USB memory 368, and the USB memory 368 is inserted into the USB port 416 and the program is transferred to the SSD 408. Alternatively, the program may be transmitted to the computer 360 via the network 370 and stored in the SSD 408.

[0083] The program is loaded into RAM 406 when executed. Of course, the source program may be input using keyboard 366, monitor 362, and mouse 364, and the compiled object program may be stored in SSD 408. In the case of a program that runs on a virtual machine, a program that functions as a virtual machine must be installed in advance on computer 360.

[0084] The CPU 400 reads a program from the RAM 406 according to an address indicated by an internal register called a program counter (not shown), interprets the instructions, reads data required to execute the instructions from the RAM 406, the SSD 408, or another device according to the address specified by the instruction, and executes the processing specified by the instruction. The CPU 400 stores the execution result data at an address specified by the program, such as the RAM 406, the SSD 408, or a register within the CPU 400. At this time, the value of the program counter is also updated by the program. The computer program may be loaded directly into the RAM 406 from the USB memory 368 or via a network. Note that some tasks (mainly numerical calculations) of the program executed by the CPU 400 may be executed by the GPU 402 according to instructions included in the program or according to analysis results obtained when the CPU 400 executes the instructions, thereby enabling the quantum Monte Carlo calculation method described above to quickly output results of combinatorial optimization problems.

[0085] A program that cooperates with computer 360 to implement each step of the quantum Monte Carlo calculation method, each function of the quantum Monte Carlo calculation device, and the method for solving a path optimization problem includes a plurality of instructions written and arranged to cause computer 360 to operate to implement those functions. Some of the basic functions required to execute the instructions are provided by an operating system (OS) or third-party programs running on computer 360, or by modules of various toolkits installed on computer 360. Therefore, the program does not necessarily include all of the functions required to implement the system and method of this embodiment. The program need only include instructions for executing the operations of the intended devices and their components by statically linking appropriate functions or program routines in a controlled manner to achieve the desired results, or by dynamically linking to those functions during program execution. The method for operating computer 360 in this manner is well known. Therefore, a description of the method for operation will not be repeated here.

[0086] 2. Second Embodiment 2.1. Configuration In the first embodiment, as shown in FIG. 5, both the temperature and the transverse magnetic field are changed to ⊥is arbitrarily varied. However, this disclosure is not limited to such an embodiment. This second embodiment is an example in which the temperature T is constant (line 450), as shown in FIG. 11 . Quantum computers currently being commercialized operate only at extremely low temperatures. Therefore, it is difficult to manipulate the temperature. On the other hand, it is relatively easy to control the transverse magnetic field. Therefore, this second embodiment is considered to be executable in quantum computers. In this embodiment, the transverse magnetic field is controlled according to a function that rapidly decreases from a constant value ("1") as the Monte Carlo steps progress, as shown by curve 454, and then converges to 0. Various functions with such characteristics are conceivable. In this embodiment, a function with the same form as 1 / 100 (1.01-x) is used as this function (the Monte Carlo steps are normalized to the range of 0 to 1). FIG. 11 shows the J corresponding to the Monte Carlo steps. ⊥ The change in is also shown (curve 452).

[0087] This second embodiment also uses the same basic concept as the first embodiment. However, the configuration of this second embodiment differs in some details from the first embodiment. Figure 12 shows, in flowchart form, the control structure of the program executed in this second embodiment. This program realizes processing equivalent to step 154 ​​in Figure 3 in the first embodiment. In other words, this program corresponds to the program shown in Figure 4 in the first embodiment, but has a slightly different configuration from that shown in Figure 4.

[0088] 12, program 470 differs from the program shown in Fig. 4 in that it includes step 480 instead of step 206 in Fig. 4. In step 480, step 482 is repeatedly executed a predetermined number of times.

[0089] Step 480 includes step 500 of repeating step 502 for the number of Trotter layers, and step 504 of updating (decreasing) the value of the transverse magnetic field (Γ) after completing step 500. In step 504, the temperature (T) is not updated but is kept constant.

[0090] Step 502 includes steps 250, 252, and 254, which are similar to step 232 in Fig. 4. However, step 502 differs from step 232 in that step 502 includes step 520 instead of step 256 in Fig. 4. Details of step 520 are shown in Fig. 13.

[0091] 13, the program of step 520 includes step 550, which branches the control flow depending on whether the energy E1 calculated in step 252 is smaller than the minimum energy value Emin found in the processing up to that point. Note that "smaller or not" in step 550 means "less than or equal to." However, an embodiment is also possible in which the determination in step 550 is "equal to or less than or equal to."

[0092] Step 520 further assigns the value of the new energy E1 to the minimum energy value Emin when the determination in step 550 is affirmative, and calculates the spin arrangement |α ｍｉｎ 〉 into the current spin configuration |α Ｋ >, and following step 552, the current spin configuration |α Ｋ 〉' as a state candidate |α Ｋ and terminating execution of the program.

[0093] This program further includes step 310 for sampling a random number ξ in the interval [0, 1] when the determination in step 550 is negative, and setting the value of the random number ξ sampled in step 310 to e （－λΔＥＴ） and step 554, which branches the flow of control depending on whether or not the exponent "ΔE" is smaller than the exponent "ΔE1" in step 312 of the first embodiment in FIG. 4. If the determination in step 554 is affirmative, control proceeds to step 556 described above. If the determination in step 554 is negative, execution of this program ends.

[0094] The delivery plan optimization system using the program according to the second embodiment operates in almost the same manner as the system according to the first embodiment. The differences between the two systems are based on the structural differences described above. Therefore, the differences in operation based on these differences will not be repeated here.

[0095] 2.2 Benchmarks With reference to FIGS. 14 to 17, the effects of the quantum Monte Carlo calculations in the second embodiment will be described based on the results of the benchmark problems described above.

[0096] FIG. 14 shows the results (curve 580) of the second embodiment for the benchmark problem data set G9, along with the results (curve 582) of the first embodiment. This processing was performed on data set G9 with the temperature T fixed at 0.004 and the initial value of the transverse magnetic field Γ set to 1.0. FIG. 14 also shows an enlarged graph 584 of a portion of curve 580 and curve 582. Enlarged graph 584 shows the portion of curve 580 and curve 582 that transitions from a steep rise to a gentler waveform.

[0097] Similarly, for data set G20, the results of the second embodiment (curve 600) are shown in Figure 15 together with the results of the first embodiment (curve 602). This processing was performed on data set G20 with the temperature T fixed at 0.002 and the initial value of the transverse magnetic field Γ set to 1.0. In Figure 15, the portions of curves 600 and 602 that change from a rising edge to a smooth waveform are shown as enlarged graphs 604.

[0098] Comparing curves 580 and 582 in FIG. 14 and curves 600 and 602 in FIG. 15 reveals that the final convergence value in the second embodiment is substantially the same as that in the first embodiment. Meanwhile, in the example shown in FIG. 14 , it can be seen that the rise of curve 580 is significantly faster than that of curve 582, and in the example shown in FIG. 15 , it can be seen that the rise of curve 600 is even faster than that of curve 602. That is, in the second embodiment, a value close to the final convergence value is already obtained at a fairly early stage of the Monte Carlo step. Therefore, according to the second embodiment, the convergence speed is significantly improved. As a result, the second embodiment can be used when a rough solution needs to be known quickly.

[0099] The improvement in convergence speed becomes more apparent as the number of nodes in the graph to be processed increases. Figure 16 shows the results of quantum Monte Carlo calculations according to the second embodiment (curve 620) and the results of the conventional SA method (curve 622) for dataset G34. Figure 16 also shows a logarithmic plot 624, in which the horizontal axis represents the logarithmic representation of the number of Monte Carlo steps for curves 620 and 622. In the logarithmic plot 624, curve 630 corresponds to curve 620. Curve 632 corresponds to curve 622.

[0100] 16, the difference in convergence speed between the method of the second embodiment and the conventional SA method is clear. In particular, looking at curve 632 of logarithmic plot graph 624, a clear stagnation is observed during the calculation in the conventional SA method. However, no such stagnation is observed in curve 630 corresponding to the second embodiment.

[0101] Similarly, Figure 17 compares the results (curve 650) of the second embodiment with the results (curve 652) of the first embodiment for data set G34. Figure 17 also shows a logarithmic plot 654, in which the horizontal axis represents the logarithmic representation of the number of Monte Carlo steps for curves 650 and 652. In logarithmic plot 654, curve 660 corresponds to curve 650. Curve 662 corresponds to curve 652.

[0102] 17 shows that the convergence speed of the quantum Monte Carlo calculation according to the second embodiment is faster than that of the calculation according to the first embodiment. The following should be noted in particular. That is, in the calculation according to the first embodiment, as shown by curve 662 of the logarithmic plot graph 654, stagnation is observed midway, although not as severe as in the SA method. In contrast, in the calculation according to the second embodiment, as shown by curve 660, no stagnation is observed midway. Therefore, it can be seen that the quantum Monte Carlo calculation method according to the second embodiment has a faster convergence speed than the first embodiment.

[0103] As described above, according to each embodiment of the present disclosure, a quantum Monte Carlo calculation method that can shorten the calculation time compared to conventional methods can be realized. As a result, according to this disclosure, a high-speed approximate solution method for problems that require a huge amount of calculation, such as a route optimization problem, can be provided.

[0104] Each process (each function) in each of the above-described embodiments is realized by a processing circuit (circuitry) including one or more processors. The processing circuit may be configured by an integrated circuit that combines one or more memories, various analog circuits, and various digital circuits in addition to the one or more processors. The one or more memories store programs (instructions) that cause the one or more processors to execute each of the processes. The one or more processors may execute each of the processes according to the program read from the one or more memories, or according to a logic circuit designed in advance to execute each of the processes. The processor may be any of various processors suitable for computer control, such as a CPU, GPU, DSP (Digital Signal Processor), FPGA (Field-Programmable Gate Array), ASIC (Application Specific Integrated Circuit), and DPU (Data Processing Unit). The physically separated processors may cooperate with each other to execute the processes. For example, the processors installed in the physically separated computers may cooperate with each other via a network such as a LAN (Local Area Network), a WAN (Wide Area Network), or the Internet to execute the processes. The program may be installed into the memory from an external server device or the like via the network, or may be distributed in a state stored on a recording medium such as a CD-ROM (Compact Disc Read-Only Memory), a DVD-ROM (Digital Versatile Disc Read-Only Memory), or a semiconductor memory, and then installed into the memory from the recording medium.

[0105] The above-described embodiment relates to a two-dimensional Ising model. However, this disclosure is not limited to such an embodiment. It can also be applied to a three-dimensional Ising model, and in principle, it can be applied to any calculation method that performs Monte Carlo calculations using the Suzuki-Trotter decomposition.

[0106] The embodiments disclosed herein should be considered in all respects as illustrative and not restrictive. The scope of the disclosure is not defined by the detailed description of the disclosure, but by the claims of the appended claims, and is intended to include all modifications within the meaning and scope of the claims.

[0107] 50, 90 Ising model 60, 62, 64 Node 70, 72 Spin 100, 102, 104 Trotter layer 280, 450 Line 282, 330, 332, 452, 454, 580, 582, 600, 602, 620, 622, 630, 632, 650, 652, 660, 662 Curve 350 Computer system 360 Computer 362 Monitor 364 Mouse 366 Keyboard 368 USB memory 370 Network 400 CPU 402 GPU 404 ROM 406 RAM 408 SSD 410 Bus 412 Display control unit 414 Input / output I / F 416 USB port 418 Network I / F 470 Program 584, 604 Enlarged graph 624, 654 Logarithmic plot graph

Claims

1. A quantum Monte Carlo calculation method for calculating a minimum energy state that gives a minimum energy, which is a minimum value of energy in an Ising model, comprising: a step in which a computer prepares a predetermined number of spin arrays for the Ising model; a step in which a computer sets an initial state of the minimum energy state for the Ising model and an initial value of the minimum energy; and a Monte Carlo calculation step in which a computer repeats a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until a termination condition is met, wherein the Monte Carlo calculation step includes: for each of the predetermined number of spin arrays, a step in which a computer generates a new state candidate in which the spin state of the spin array to be processed has changed; and a step in which a computer performs an update process for the state of the spin array to be processed and the minimum energy state based on the energy for the state candidate, the energy for the spin array to be processed, and the minimum energy obtained up to a current time.

2. The energy of the Ising model is composed of a sum of a first energy due to interactions within the spin array to be processed and a second energy due to interactions between the spin array to be processed and the spin array that gives the minimum energy obtained up to the present time, and the step of performing the update process includes a first determination step in which a computer determines whether or not a difference between the first energy related to the state candidate and the first energy related to the spin array to be processed is smaller than 0; and a first selective processing step in which, in response to a positive determination in the first determination step, the computer selectively executes a process of updating the state of the spin array to be processed with the state candidate based on the first energy related to the state candidate and the minimum energy obtained up to the present time, and updating the minimum energy state and the minimum energy with the state candidate and the first energy related to the state candidate, respectively, and a process of updating the state of the spin array to be processed with the state candidate; 2. The quantum Monte Carlo calculation method according to claim 1, further comprising: a second selective processing step in which, in response to a negative determination in the first determination step, a computer selectively executes a process of updating or not updating the state of the spin array to be processed with the state candidate, depending on whether a difference between the energy related to the state candidate and the energy related to the spin array to be processed is smaller than 0.

3. The quantum Monte Carlo computational method of claim 2, wherein a contribution of the first energy in the sum of the first energy and the second energy is greater than a contribution of the second energy.

4. The quantum Monte Carlo computational method of claim 3, wherein said contribution of said first energy is one and said contribution of said second energy is zero.

5. A quantum Monte Carlo calculation method according to any one of claims 2 to 4, wherein the first selective processing step includes: a second determination step in which, in response to a positive determination in the first determination step, a computer determines whether or not the first energy related to the state candidate is smaller than the first energy related to the spin array to be processed; a first updating step in which, in response to a positive determination in the second determination step, a computer updates the state of the spin array to be processed with the state candidate and updates the minimum energy state and the minimum energy with the state candidate and the first energy related to the state candidate, respectively; and a second updating step in which, in response to a negative determination in the second determination step, a computer updates the state of the spin array to be processed with the state candidate.

6. A quantum Monte Carlo calculation method as described in claim 5, wherein the second selective processing step includes: a third determination step in which, in response to a negative determination in the first determination step, the computer determines whether or not a difference between the energy associated with the state candidate and the energy associated with the spin array to be processed is equal to or greater than 0; and a third update step in which, in response to a positive determination in the third determination step, the computer performs processing to update the state of the spin array to be processed with the state candidate.

7. A quantum Monte Carlo calculation method as described in claim 6, wherein the third updating step includes a step in which a computer updates the state of the spin array to be processed with the state candidate according to a probability determined by a difference between the first energy related to the state candidate and the first energy related to the spin array to be processed.

8. The quantum Monte Carlo calculation method according to any one of claims 1 to 4, wherein the Monte Carlo calculation step is executed in parallel for the predetermined number of spin arrays.

9. A quantum Monte Carlo calculation method according to any one of claims 1 to 4, wherein the step of generating the new state candidate includes a step in which a computer randomly selects one of the nodes of the spin array to be processed and inverts the spin of the node.

10. The quantum Monte Carlo calculation method according to claim 1, wherein the energy of the Ising model is composed of a sum of a first energy due to interactions within the spin array to be processed and a second energy due to interactions between the spin array to be processed and the spin array that gives the minimum energy obtained up to a current time, and the step of performing the updating process includes an updating step in which a computer executes a process of updating the minimum energy with the energy of the spin array to be processed in response to the first energy for the state candidate being smaller than the minimum energy obtained up to a current time, a process of updating the state of the spin array to be processed with the state candidate, and a process of updating the spin array with the state candidate in accordance with a probability determined by a difference between the energy for the state candidate and the energy for the spin array to be processed.

11. A quantum Monte Carlo calculation device for calculating a minimum energy state that gives a minimum energy in an Ising model, comprising: a Trotter layer preparation unit that prepares a predetermined number of spin arrays for the Ising model; an initial state setting unit that sets an initial state of the minimum energy state for the Ising model and an initial value of the minimum energy; and a Monte Carlo calculation unit that repeats a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until a termination condition is met, wherein the Monte Carlo calculation unit includes: a state candidate generation unit that generates, for each of the predetermined number of spin arrays, a new state candidate in which the spin state of the spin array to be processed has changed; and an update processing unit that performs an update process of the state of the spin array to be processed and the minimum energy state, for each of the predetermined number of spin arrays, based on the energy for the state candidate, the energy for the spin array to be processed, and the minimum energy obtained up to a current time.

12. A computer program causing a computer to execute a quantum Monte Carlo calculation method for calculating a minimum energy state that gives a minimum energy in an Ising model, the quantum Monte Carlo calculation method including: a step in which a computer prepares a predetermined number of spin arrays related to the Ising model; a step in which a computer sets an initial state of the minimum energy state related to the Ising model and an initial value of the minimum energy; and a Monte Carlo calculation step in which a computer repeats a Monte Carlo calculation process for searching for the minimum energy state of the Ising model until a termination condition is met, the Monte Carlo calculation step including: for each of the predetermined number of spin arrays, a step in which a computer generates a new state candidate in which the spin state of the spin array to be processed has changed; and a step in which a computer performs an update process for the state of the spin array to be processed and the minimum energy state based on the energy related to the state candidate, the energy related to the spin array to be processed, and the minimum energy obtained up to the current time.

13. A method for solving a route optimization problem, comprising: a step of converting, by a computer, a route optimization problem formulated in a predetermined format into an Ising model; a step of determining, by a computer, a solution for the Ising model by applying the quantum Monte Carlo calculation method described in claim 1 to a spin arrangement derived from the Ising model; and a step of converting, by a computer, the solution into a solution for the route optimization problem.

Images