COMPUTING APPARATUS AND COMPUTING METHOD

Abstract

An object of the invention is to provide a computing technology which can operate at room temperature and have a sufficient performance for combinatorial optimization problems that need an exhaustive search. In a local-field response method in which spins being variables respond to local effective magnetic fields, a time axis is discretely treated. When spins respond to effective magnetic fields, the effective magnetic fields are determined sequentially from the site having the small magnitude of a spin, and spins respond to the fields in order. When the sign of a spin is inverted, the information is reflected in the subsequent process of determining the effective magnetic fields for other sites. Thus, a many-body effect due to quantum entanglement is phenomenologically incorporated.

Description

CROSS-REFERENCE TO RELATED APPLICATION

The present application claims priority from Japanese application JP 2017-204990, filed on Oct. 24, 2017, the content of which is hereby incorporated by reference into this application.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to a computing technology which is able to perform computation at a high speed with respect to combinatorial optimization problems that need an exhaustive search.

2. Related Art

As being representative as the word “IoT” (Internet of Things), various things are connected to the Internet in the present days; information is collected from the things; and the things are controlled by the collected information. In the control, an optimal solution is found from among many choices and is performed. Extremely speaking, the information technology in the present days can be said to search an optimal solution.

In this situation, a quantum annealing, or adiabatic quantum computing, has recently received attention. In this method, a problem is set such that the ground state of a certain physical system is a solution and the solution is obtained through finding the ground state. Let the Hamiltonian of a physical system in which a problem is set be Ĥ_p. At the beginning of computation, however, Hamiltonian is not Ĥ_p but Ĥ₀ for which the ground state is easily prepared. Next, the Hamiltonian is transformed from Ĥ₀ to Ĥ_p with a sufficient period of time. When the transformation takes a sufficient period of time, the system continues to stay in the ground state, and the ground state (solution state) for Hamiltonian Ĥ_p is finally obtained. This is a principle of the quantum annealing.

A ground state-searching method in which a physical system called Ising spin glass is used can be applicable even to a problem called NP-hard. Meanwhile, a highly difficult problem in combinatorial optimization problems belongs to the NP-hard. Moreover, problems classified into “P” and problems classified into “NP” in the computational complexity theory all can be transformed to an NP-hard problem. Therefore, if the quantum annealing is applied to the Ising spin glass system, almost all the combinatorial optimization problems can be solved, and the most important issue of the information technology is solved.

Another reason why the quantum annealing receives attention is robustness with respect to decoherence. In a quantum computer, quantum coherence must be preserved over a computing time. However, in the quantum annealing, the condition is relaxed. If the ground state is maintained, a correct solution is obtained. The quantum coherence is not necessarily to be maintained. In a current technology level, it is difficult to establish a pure quantum system. Therefore, the quantum coherence is hardly kept over the computing time. For this reason, the quantum annealing has received attention. However, there is a drawback even in the quantum annealing. The quantum annealing can only be realized in a superconducting magnetic flux qubit system as it stands, and thus, a cryogenic cooling apparatus is needed. The need for an extremely low temperature is an issue of achieving a practical computer.

To solve the issue, the local-field response method was proposed as described in the following (PTLs 1-3 and NPL 1). First, let us review the quantum annealing again. The concept of the annealing lies regardless of quantum or classical by its nature. The quantum annealing is devised to improve the performance of the classic annealing using quantum properties. The reason why the quantum annealing does not require the quantum coherence over the computing time and only requires the ground state being maintained comes from the wide concept of the annealing. As just described, the concept of the annealing is so wide that another method to use quantum properties may exist, different from the quantum annealing.

The local-field response method has been disclosed from such a point of view. In this method, similarly to the quantum annealing, a transverse field is applied to a spin system at time t=t₀, which is a computing device, and the magnetic field is gradually decreased to obtain a solution at time t=τ. The computing device itself is classical, and quantum-mechanical information is added to the response of spins to the magnetic field. The method works on a classical machine at room temperature, and therefore, it can solve the issue of cryogenic cooling in the quantum annealing. In PTLs 1-3 and NPL 1, quantum effects were introduced in the response function, where the response function which has an average quantum effect was determined empirically or based on the results solved for similar problems. In NPL 2, solution accuracy was improved compared to the methods in PTLs 1-3 and NPL 1 by phenomenologically incorporating the properties of the linear superposition in quantum mechanics. However, the property of quantum entanglement which is another important property in quantum mechanics has not been sufficiently incorporated.

CITATION LIST

Patent Literature

• PTL 1: WO 2015/118639
• PTL 2: WO 2016/157333
• PTL 3: WO 2016/194221

Non-Patent Literature

• NPL 1: T. Tomaru, “Quasi-Adiabatic Quantum Computing Treated with c-Numbers Using the Local-Field Response,” J. Phys. Soc. Jpn. 85, 034802 (2016).
• NPL 2: T. Tomaru, “Improved Local-field response Method Working as Quasi-Quantum Annealing,” J. Phys. Soc. Jpn. 86, 054801 (2017).

SUMMARY OF INVENTION

Technical Problem

As described above, the quantum annealing requires a cryogenic cooling apparatus to use a superconducting magnetic flux qubit system. Meanwhile, the local-field response method operates at room temperature, but quantum entanglement which is an important property in quantum mechanics is not sufficiently incorporated, which limits the performance. Thus, an object of the invention is to provide a computing apparatus and a computing program which can operate at room temperature and have a sufficient performance for a difficult assignment such as an issue that needs an exhaustive search.

Solution to Problem

In a local-field response method in which spins being variables respond to local effective magnetic fields, a time axis is discretely treated; When spins respond to effective magnetic fields, the effective magnetic fields are determined sequentially from the site having the small magnitude of a spin, and the spins respond to the fields in order. When the sign of the spin is inverted, the information is reflected in the subsequent process of determining the effective magnetic fields for other sites. Thus, a many-body effect due to quantum entanglement is phenomenologically incorporated. More specific descriptions are provided in the following.

A specific figure is a computing apparatus which includes a computing unit, a storage unit, and a control unit, and performs computation under the control of the control unit while transferring data between the storage unit and the computing unit, or a computing method using the computing apparatus. N variables s_j^z (j=1, 2, . . . , N) take a range of −1≤s_j^z≤1, and an assignment is set with coefficients g_j indicating local terms and coefficients J_kj (k, j=1, 2, . . . , N) indicating inter-variable interactions. Time is divided into m, and the computing unit discretely performs computation from t=t₀ (t₀=0) to t_m (t_m≤τ). Herein, N and m are natural numbers.

Variables B_eff,j^z(t_i) and s_j^z(t_i) at each time t_i (i=1, 2, . . . , m) are determined in this order. B_eff,j^z(t_i) is a function of s_k^z(t_i−1), J_kj, g_j, and t_i. s_j^z(t_i) is a function of B_eff,j^z(t_i) and t_i. Initial values at time t₀ are set as B_j^z(t₀)=0 and s_j^z(t₀)=0. For determining B_eff,j^z(t_i) and s_j^z(t_i) at time t_i (i=1, 2, . . . , m), first, s_j^z(t_i−1) are put in descending order such that |s_m1^z(t_i−1)|≤|s_m2^z(t_i−1)|≤s_m3^z(t_i−1)|≤ . . . ≤|s_mN^z(t_i−1)|. Then, B_eff,m1^z(t_i) and s_m1^z(t_i) at site m₁ are determined at the first time, and s_m1^z(t_i−1) is set to be s_m1^z(t_i−1)=sgn(s_m1^z(t_i))|s_m1^z(t_i−1)|. Next, B_eff,m2^z(t_i) and s_m2^z(t_i) at site m₂ are determined and s_m2^z(t_i−1) is set to be s_m2^z(t_i−1)=sgn(s_m2^z(t_i))|s_m2^z(t_i−1)|. The computation at site m₃ is performed similarly, and the computation up to site m_N is performed similarly for the computation at time t_i. As the time step progresses from t=t₀ to t=t_m, the variable s_j^z approaches −1 or 1, and a solution is determined as s_j^zfd=−1 if s_j^z<0 and as s_j^zfd=1 if s_j^z>0.

As described, the same computation has been performed up to the site m_N for the computation at time t_i, but there is an option that the same computation is performed up to m_x (herein, 0<x≤N), and that after site m_x all remaining sites are processed in an individual and parallel manner similarly to the original local-field response method for the computation at time t_i. Here, x is a natural number.

Advantageous Effects of Invention

The present method operates on a classical machine. Extremely low temperature is not needed. In addition, there is no need to take quantum coherence into consideration. As a result, usable resources are expanded, and electric circuits can also be used. Moreover, solution accuracy is improved and a computation time is shortened through phenomenologically incorporating the effect of quantum entanglement. Thanks to these properties, it is possible to realize a practical computing apparatus which can solve difficult problems with high solution accuracy.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a diagram schematically illustrating a principle of an embodiment;

FIG. 2 is a flowchart illustrating an example of an algorithm in accordance with the first embodiment;

FIG. 3 is a graph plotting specific values of the response function r_b;

FIG. 4 is a flowchart illustrating an example of an algorithm in accordance with the second embodiment;

FIG. 5 is a flowchart illustrating an example of an algorithm in accordance with the third embodiment;

FIG. 6 is a flowchart illustrating another example of the algorithm in accordance with the third embodiment;

FIG. 7 is a flowchart illustrating an example of an algorithm in accordance with the fourth embodiment;

FIG. 8 is a flowchart illustrating an example of an algorithm in accordance with the fifth embodiment;

FIG. 9 is a diagram illustrating an example of an algorithm in accordance with the sixth embodiment; and

FIG. 10 is a block diagram illustrating an example of a configuration of a computing apparatus in accordance with the seventh embodiment.

DESCRIPTION OF EMBODIMENTS

Embodiments will be described in detail using the drawings. However, the content of the embodiments described below is not interpreted in a way of limiting the invention. A person skilled in the art can easily understand that the specific configuration may vary in a scope not departing from the idea and the spirit of the invention.

Portions having the same or similar functions in the configuration of the invention described below will be represented with the same symbol and commonly used in different drawings, and the redundant description will be omitted.

In a case where there are elements having similar or corresponding functions, the elements may be attached to the same symbol with different suffixes. However, in a case where there is no need to distinguish plural elements, the suffixes may be omitted.

First Embodiment

In the first embodiment, we will give the basic principle, starting from a quantum-mechanical description and transforming it into a classical form.

FIG. 1 schematically show the principle of the embodiment. A basic framework is the same as the local-field response method disclosed in PTL 1 and NPL 1. A transverse field is applied at t=0 to make spins directed in one direction. Thereafter, the transverse field is gradually decreased, and Hamiltonian is set to the problem Hamiltonian at t=τ. Spins time-evolve in response to the local effective magnetic field which is applied to each spin at each time.

Let the problem Hamiltonian and the Hamiltonian at t=0 be Eqs. (1) and (2), respectively.

$\begin{array}{cc}{\hat{H}}_p=-\sum _{i>j}J_{\mathrm{ij}}{\hat{\sigma }}_i^z{\hat{\sigma }}_j^z-\sum _jg_j{\hat{\sigma }}_j^z& \left(1\right)\\ {\hat{H}}_0=-\gamma \sum _j{\hat{\sigma }}_j^z& \left(2\right)\end{array}$

Let the Hamiltonian at time t be Eq. (3).

$\begin{array}{cc}\hat{H}\left(t\right)=\left(1-\frac{t}{\tau }\right){\hat{H}}_0+\frac{t}{\tau }{\hat{H}}_p& \left(3\right)\end{array}$

Herein, τ is a computation time. From an analogy with a one-spin system, the effective magnetic field at site j is given by B̂_eff,j=−∂Ĥ/∂σ̂_j.

$\begin{array}{cc}{\hat{B}}_{\mathrm{eff},j}\left(t\right)=\left(\left(1-\frac{t}{\tau }\right)\gamma ,0,\frac{t}{\tau }\left(\sum _{k\left(\ne j\right)}J_{\mathrm{kj}}{\hat{\sigma }}_k^z+g_j\right)\right)& \left(4\right)\end{array}$

The local-field response method of this embodiment operates on a classical machine in which expectation value <σ̂_j> is regarded as a spin variable. As seen in Eqs. (1) and (2), <σ̂_j> and <B̂_eff,j> consists of only x and z components. Therefore, a response function r_b(t) is defined only by the x and z components as described by Eq. (5).

{circumflex over (σ)}_j^z(t)/{circumflex over (σ)}_j^x(t)=r_b(t)·{circumflex over (B)}_eff,j^z(t)/{circumflex over (B)}_eff,j^x(t)  (5)

A spin direction is determined based on the response function. When the spin system is classical, the response of each spin is determined only by the effective magnetic field at each site, and the response function is r_b(t)=1.

However, there is a non-local correlation (quantum entanglement) in quantum mechanics, and generally r_b(t)≠1. As mentioned already, Eq. (5) has been transformed to a classic form by taking the expectation value, but a quantum effect is incorporated in r_b(t)≠1. A value of r_b(t) is determined empirically or through quantum-mechanical prior calculations for similar problems. The quantum effect herein is an average one because it is empirically determined or it is based on the results for similar problems. The local-field response method can incorporate quantum effects through various methods. The method through r_b(t)≠1 is an example. Note that setting r_b(t)=1 without the quantum effect is also allowed in the local-field response method. The local-field response method itself can operate even when the quantum effect is not incorporated.

Four variables <σ̂_j^z(t_i)>, <σ̂_j^x(t_i)>, <B̂_eff,j^z(t_i)>, and <B̂_eff,j^x(t_i)> in Eq. (5) are expectation values, and therefore, classical quantities. For this reason, let us change the quantum-mechanical description to a classical-physics one, i.e., <σ̂_j^x(t_i)>→s_j^x(t_i), <σ̂_j^z(t_i)>→s_j^z(t_i), <B̂_eff,j^x(t_i)>→B_eff,j^x(t_i), <B̂_eff,j^z(t_i)>→B_eff,j^z(t_i). With the description change, Eq. (5) is modified to Eq. (6).

(s_j^z(t)/s_j^x(t))=r_b(t)·(B_eff,j^z(t)/B_eff,j^x(t))  (6)

The time-evolution in the local-field response method is discrete. B_eff,j^z(t_i) at time t_i is determined from s_k^z(t_i−1) at time t_i−1 in accordance with Eq. (4). s_j^z(t_i) at time t_i is determined from B_eff,j^z(t_i) at time t_i in accordance with Eq. (6). This procedure is repeated. FIG. 2 summarizes the procedure as a flowchart.

Step 101 in FIG. 2 indicates a start point of the algorithm and sets initial values. Step 102a determines B_eff,j^z(t_i) using s_k^z(t_i−1) at time t_i−1, and determines the transverse field intensity B_eff,j^x(t_i) depending on time. Step 103 determines tan θ=r_b(t)·B_eff,j^z(t_i)/B_eff,j^x(t_i) corresponding to the spin direction using B_eff,j^z(t_i)/B_eff,j^x(t_i) and the response function r_b(t), and it determines s_j^z(t_i)=r_s(t)·sin θ using a parameter r_s(t) representing the magnitude of the spin. r_s(t) is defined by r_s(t)²=s_j^x(t_i)²=s_j^z(t_i)², and it is determined empirically or through prior calculations similarly to r_b(t). A specific example will be described in the following second paragraph. Steps 102a and 103 in FIG. 2 are a set of repeated computation. The set is repeated until t=t_m (≤τ). At t=t_m, a solution s_j^zfd is determined as s_j^zfd=1 if s_j^z(t)>0 and s_j^zfd=−1 if s_j^z(t_m)<0 (Steps 201 and 202). Herein, the reason why t_m≤τ is assumed is that many cases obtain converged solutions before t=τ in the time-evolution.

Step 103 can be generally written using a function f such as s_j^z(t_i)=f(B_eff,j^z(t_i), t_i), and f(B_eff,j^z(t_i), t_i)=r_s(t)·sin{arctan(r_b(t)·B_eff,j^z(t_k)/B_eff,j^x(t_k))}. r_s(t) is a parameter expressing the magnitude of the spin and satisfies 0≤r_s(t)≤1. Thus, −1≤s_j^z(t_i)≤1. Here, when r_b(t)=1 and r_s(t)=1, then the system is purely classical.

FIG. 3 shows an example of prior computations for the response function r_b(t). The example shows a case where J_ij and g_j are determined by uniform random numbers of [−5, 5] in an 8-bit system. The figure shows the results of 100 problems and contains 800 points (100 problems×8 bits). Here, B_zx≡<B̂_eff,j^z>/<B̂_eff,j^x> and s_zx≡<σ̂_j^z>/<σ̂_j^x>. The points indicate values computed exactly quantum-mechanically. The response functions largely scatter owing to the non-local correlation in quantum mechanics. Open circles indicate average values, where the horizontal axis is divided into 40 regions and the points in each region are averaged. The averaged response function has smooth dependence on B_zx. Because being smooth, the response function can be expressed with several parameters. NPL 1 discloses a technique of expressing the smooth response function using four parameters. A solid line r_b⁰(t) in FIG. 3 is obtained with the technique. Another parameter r_s(t) is also obtained from the same four parameters.

Hitherto, an example of the response function has been described in FIG. 3, and a method of determining r_b⁰(t) and r_s(t) has been described with reference to NPL 1. However, the method of determining the response functions r_b⁰(t) and r_s(t) is not limited to the above description, and various methods may be employed and the response functions can also be empirically determined. In the following embodiment, r_b⁰(t) is not limited to the solid line in FIG. 3, but assumed to be a response function in which quantum effects are incorporated on average.

Here, similarly to r_b(t), r_s(t) is also possible to set r_s(t)=1. While the range of r_b(t) is −∞<r_b(t)<∞, r_s(t) is 0≤r_s(t)≤1. Therefore, the influence of assuming r_s(t)=1 on a final solution is smaller than that of assuming r_b(t)=1. Thus, it is effective to set r_s(t)=1 when prior information for determining r_s(t) is insufficient.

Second Embodiment

The first embodiment has been described that quantum effects can be averagely incorporated through r_b(t)≠1. However, quantum effects depend on problems and vary with time. Quantum effects cannot be sufficiently incorporated only through averaged quantities. This embodiment describes a method of phenomenologically incorporating a quantum entanglement related-quantum effect in the formulation depending on the spin state at each time.

The influence of quantum entanglement appears as a many-body effect. When quantum entanglement is large, if a certain spin is inverted (its sign is inverted), another spin is simultaneously inverted with high probability. In the algorithm of FIG. 2, the effective magnetic field B_eff,j^z(t_i) at t=t_i is calculated site by site using s_j^z(t_i−1) at t=t_i−1. The calculation is performed independently site by site and it is in a one-body approximation. Therefore, simultaneous inversion of spins has not been sufficiently taken into consideration. For this reason, we take simultaneous inversion of spins into consideration in accordance with the algorithm in FIG. 4.

In the vicinity of the spin inversion, s_j^z≈0. Therefore, the spin that has a smaller value of |s_j^z| has higher probability of inversion. Thus, let us first put |s_j^z(t_i−1)| at each time t_i−1 in descending order.

Let m₁, m₂, m₃, . . . be the site number of |s_j^z(t_i−1)| in the descendent order as illustrated in Step 111 in FIG. 4. B_eff,m1^z(t_i) at site m₁ is calculated using Eq. (4) (Step 102a), and s_m1^z(t_i) is calculated using Eq. (6) (Step 103). If there is no sign inversion in the time evolution from s_m1^z(t_i−1) to s_m1^z(t_i), the next procedure is to compute B_eff,m2^z(t_i) and s_m2^z(t_i) at site m₂. If the sign is inverted, s_m1^z(t_i−1) is changed to −s_m1^z(t_i−1), and after that, B_eff,m2^z(t_i) and s_m2^z(t_i) are computed. In other words, the sign of s_m1^z(t_i−1) is changed into the sign of s_m1^z(t_i) as s_m1^z(t_i−1)→sgn(s_m1^z(t_i))·|s_m1^z(t_1-1)| (Step 112). Here, sgn denotes a sign function which returns 1 or −1 in accordance with the sign of the argument.

According to this procedure, the change for site m₁ at time t_i is immediately reflected on the time evolution for site m₂. When sites m₁ and m₂ are entangled, if the spin at site m₁ is inverted, the spin at site m₂ is inverted with high probability. Such an effect is incorporated as the change of the sign described herein. In other words, the effect of quantum entanglement is incorporated phenomenologically.

The computed result of s_m2^z(t_i) is also processed similarly to s_m1^z(t_i) i.e., s_m2^z(t_i−1)→sgn(s_m2^z(t_i))·|s_m2^z(t_i−1)| (Step 112). By repeating the similar process for site m₃ and the following sites, we obtain N components of s_j^z(t_i) (Step 113), and the process at time t_i is completed.

Once the process at time t_i is completed, the next procedure is the process at time t_i+1, and the similar processes are repeated until time t=t_m≤τ to obtain the final solution. The final solution is s_j^zfd=1 if s_j^z(t_m)>0 and s_j^zfd=−1 if s_j^z(t_m)<0 similarly to the case of the first embodiment (Steps 201 and 202).

The spin inversion is connected to a tunneling phenomenon in a viewpoint of quantum mechanics. According to this interpretation, the treatment in this embodiment is said to take multiple tunneling into consideration.

The energy value for the obtained solution is given by Eq. (7).

$\begin{array}{cc}{H}_p\left(t_i\right)=-\sum _{k>j}J_{\mathrm{kj}}s_k^{\mathrm{zfd}}\left(t_i\right)s_j^{\mathrm{zfd}}\left(t_i\right)-\sum _jg_js_j^{\mathrm{zfd}}\left(t_i\right)& \left(7\right)\end{array}$

The local-field response method operates similarly to quantum annealing. The spin system in which the ground state is prepared at t=0 time-evolves and comes to the ground state of the problem-embedded system at t=τ ideally. The convergence of solutions is high according to the method of this embodiment. The lowest energy state during the computation is the state at t=t_m with high probability (the ground state with high probability). However, in the method of the first embodiment, the convergence of the solution is low, and the lowest energy state might be the state at t<t_m. When the first embodiment is used, the energy needs to be calculated at each time using Eq. (7) and the state corresponding to the lowest energy needs to be selected as the final solution. The amount of computation for that purpose is O(N²) for the first term in Eq. (7) and O(N) for the second term (here, O is the Landau symbol). The sum of both terms are summarized to O(N²). In contrast, if the method of this embodiment is used, because the convergence of the solution is high, the energy does not need to be calculated, and the amount O(N²) of computation can be saved.

On the other hand, the method of this embodiment rearranges the spins in Step 111. The amount of the processing is estimated as follows. If each spin is simply compared with the other spins in a repeated manner, the amount is O(N²). This is the upper limit for the processing. The lower limit is estimated as follows. Let us consider that an array of spins has already put in the descending order, and that only the adjacent spins are compared and exchanged. If there is no exchanging, the amount of computation is O(N) because each spin is compared with only one side of adjacent spins. Because a huge number of exchanges is generally rare, the actual amount of computation approaches O(N). Thus, the overhead of this embodiment is about O(N), which is smaller than the overhead O(N²) of the first embodiment.

The amount of computation for the entire local-field response method is dominated by the computation of the effective magnetic field in Step 102a. Because every N site amounts to O(N), the total amount is O(N²). Therefore, the overhead O(N) in this embodiment is negligible in a system in which N is sufficiently large.

As described above, in the example of FIG. 2, Steps 102a and 103 are processed independently and in parallel for all sites. For this reason, although the processing speed is high, quantum entanglement is not sufficiently taken into consideration. In contrast, in the example of FIG. 4 in this embodiment, the effect of quantum entanglement is phenomenologically incorporated through adding Steps 112 and 111, which makes the influence of the other spins for the relevant spin at the moment reflect at steps 102a and 103. As a result, solution accuracy is improved; the convergence of the solution is improved; and the computation time can be saved.

Third Embodiment

Quantum-mechanically, the effective magnetic field is determined based on Eq. (4). An eigenvalue of σ̂_k^z is ±1. However, because the local-field response method operates such that a spin variable s_k^z takes an expectation value <σ̂_k^z>, |s_k^z|≤1 is satisfied. For this reason, the term Σ_k(≠j)J_kjs_k^z is generally underestimated compared with g_j.

If the computation is performed while the term of Σ_k(≠j)J_kjs_k^z is underestimated, the solution accuracy is degraded. Therefore, the value of g_j is normalized with reference to the value of s_k^z. A factor c_i=(Σ_ks_k^z(t_i−1)²/N)^1/2 is multiplied to g_j to obtain g_j^norm(t_i)=c_ig_j. If g_j^norm(t_i) is set as a local term, the contributions of the terms g_j^norm(t_i) and Σ_k(≠j)J_kjs_k^z are almost equal, and the solution accuracy is improved. Here, let m (t_m≤τ) be the number of divisions in the discrete time axis, and c₁ is set as about c₁=1/m. If c₁ is simply determined in accordance with c_i=(Σ_ks_k^z(t_i−1)²/N)^1/2 and s_k^z(t₀)=0, then c₁=0. The setting of c₁=1/m is to cope with c₁=0.

FIG. 5 illustrates a flowchart containing the above processes. A difference from FIG. 4 is that Step 102a is replaced with Step 102b. In Step 102b, the process about the factor c_i=(Σ_ks_k^z(t_i−1)²/N)^1/2 is added.

FIG. 6 illustrates a modification. If a parameter c_a is introduced and the factor c_i is set as c_a·c_i, the magnitude of the factor c_i can be adjustable. FIG. 6 is such a case. As described below, the factor c_i in Step 10c in FIG. 9 is also replaced with c_a·c_i. The value of c_a is empirically determined, and is about 1 to 50.

Fourth Embodiment

In quantum-mechanical spin systems, the spins always affect other spins with each other. That is, spin σ̂_j^z at a site j affects spine σ̂_k^z at another site k, and inversely σ̂_k^z affects σ̂_j^z. Therefore, the spin σ̂_j^z affects itself through the spin σ̂_k^z at site k. That is the reason why, in quantum mechanics, the state of a spin depends not only on the state of other spins interacting with the relevant spin but also on the state of itself. A magnitude of the influence on itself through the interaction is proportional to Σ_k(≠j)J_kj². The first embodiment has described the averaged quantum effect. Because Σ_k(≠j)J_kj² are squared terms, they are left even after averaged. That is the reason why the averaged response function becomes r_b⁰(t)≠1. The open circles and the solid line in FIG. 3 show such a behavior. A detailed theory about Σ_k(≠j)J_kj² is disclosed in NPL 2.

Terms Σ_k(≠j)J_kj² that are the origin of r_b⁰(t)≠1 contain information J_kj for each problem. If the information is usable, the computation becomes more accurate. Thus, we replace the averaged response function r_b⁰(t) with r_b⁰_mod(t), where r_b⁰_mod(t) is defined by 1−r_b⁰_mod(t)=(1−r_b⁰(t))Σ_k(≠j)J_kj²/Σ_k(≠j)ave(J_kj²). Herein, ave(J_kj²) is the average of J_kj² which are used to determine r_b⁰(t). With the improvement, the response function reflects an actual problem, and the solution accuracy increases.

FIG. 7 illustrates a flowchart containing the above processes. A difference from FIG. 5 is that Step 103 is replaced with Step 103c.

Fifth Embodiment

Quantum entanglement and linear superposition are representative characteristic natures in quantum mechanics. The effect of the former quantum entanglement has been phenomenologically incorporated in up to the fourth embodiment. The effect of the latter linear superposition is phenomenologically incorporated in NPL 2. In this embodiment, the both effects are incorporated at the same time. FIG. 8 illustrates such a case.

In the embodiment illustrated in FIG. 8, the difference from FIG. 7 is that Step 102b is replaced with Step 102d.

The characteristic of linear superposition is prominent in a time region where the sign of a spin changes.

Quantum-mechanically, bands anti-cross in the vicinity of that region, and the state is linear superposition there. As a result, s_j^z(t)≈0, and correspondingly B_eff,j^z(t)≈0. In order to phenomenologically incorporate the effect, the effective magnetic fields at times t_i and t_i−1 are linearly combined to give the effective magnetic field B_eff,j^z(t_i) at t=t_i. Specifically, B_j^z0(t_i) defined by Eq. (8) is first determined using the spin values s_j^z(t_i−1) at time t_i−1.

$\begin{array}{cc}{B}_j^{z0}\left(t_i\right)=\sum _{k\left(\ne j\right)}J_{\mathrm{kj}}s_k^z\left(t_{i-1}\right)+c_ig_j& \left(8\right)\end{array}$

Next, the term at prior time is taken into consideration, and B_j^z(t_i) defined by Eq. (9) is computed.

B _j ^z(t_i)=(1−u)B_j^z0(t_i)+uB_j^z(t_i−1)  (9)

Herein, u is appropriately adjusted within 0≤u≤1 to make the solution accuracy high, typically u≈0.1. The effective magnetic field, including a transverse field and its schedule, is given by Eq. (10).

$\begin{array}{cc}{B}_{\mathrm{eff},j}\left(t_i\right)=\left(\left(1-\frac{t_i}{\tau }\right)\gamma ,0,\frac{t_i}{\tau }B_j^z\left(t_i\right)\right)& \left(10\right)\end{array}$

When the effective magnetic field is obtained, s_j^z(t) is determined using the response function r_b⁰_mod(t) in accordance with Eq. (11).

s _j ^z(t_i)/s_j^x(t_i)=r_{b mod}⁰·B_eff,j^z(t_i)/B_eff,j^x(t_i)  (11)

Sixth Embodiment

Up to the fifth embodiment, we have described the embodiments to improve the solution accuracy by adding various quantum effects, especially by adding the effect of quantum entanglement. However, the correct solution is not necessarily derived even if the above methods are used. For this reason, this embodiment will describe an auxiliary means.

Because the local-field response method is a deterministic method, the result is always the same if the same process is performed with the same initial spin values. If the initial values or the process is changed, the result may be different. Hence, we sweep the magnetic field several times while changing the initial values or the process, and we select the spin state giving the lowest energy in the process as the final solution to further improve the solution accuracy. FIG. 9 illustrates a flowchart for that case.

In FIG. 9, Processes 10a to 10d represent those described in FIG. 2 and FIGS. 4 to 8. Energy H_p^q(t_m) (q=1, 2, . . . ) is calculated in Steps 303 using spin values s_j^zfd obtained in Processes 10a to 10d, and the lowest value E_best is determined in Step 304. The spin values giving the lowest energy become the final solution.

Process 10a just performs that in FIG. 2 and FIGS. 4 to 8. In Process 10b, the sign of the result s_j^zfd obtained in Process 10a is inverted, and the sign-inverted value is set as the initial value. Strictly, the initial value is s_j^z(t₀)=s_j^zfd/div using a parameter div. The reason for dividing it with div is to sufficiently reduce the initial value, and div is set to be about m. The reason for inverting the sign is that the sign-inverted state is located farthest from the uninverted state in the spin-configuration space, and that the sign inversion gives an opportunity to escape from a local minimum if a spin state is trapped there. Process 10c sets the initial spin values as 0 similarly to Process 10a. However, a factor c_a is multiplied to a coefficient c_i such as c_a·c_i for the coefficient of the local term g_j (see the third embodiment) so as to change a balance between the interaction term and the local term. The value of c_a is empirically determined, and it is about 1 to 50. Process 10d determines the initial spin values using random numbers. Using random numbers expands possibilities of finding the solution, and the process can be simplified. Process 10d is repeatedly performed while changing the random numbers for the initial values. The solution accuracy increases as repetitions increase.

Seventh Embodiment

This embodiment is described as an algorithm, which may be executed as a software on a general computer or on a dedicated hardware. The local-field response method in the embodiments has a feature that a relatively simple computation is repeated. Therefore, it is effective that the repeated computation part is established with a dedicated hardware, and the other parts are achieved with a general-purpose device.

FIG. 10 illustrates an example of a configuration for the computation in this embodiment. FIG. 10 is similar to the configuration of a general computing apparatus, but it includes a local-field response computing device 600. The local-field response computing device 600 is the dedicated part for the computation described in the first to sixth embodiments. The other general computation is performed with a general computation device 502.

The above configuration may be constructed as a single computer. Alternatively, the configuration may be constructed from different computers connected through a network, where arbitrary parts such as a main storage device 501, a general computation device 502, a control device 503, an auxiliary storage device 504, an input device 505, and an output device 506 are connected to the network.

General computation is performed in the same procedure as in an ordinary computing apparatus. The main storage device 501 (storage unit) and the general computation device 502 (computing unit) repeatedly transfer data between them so as to progress the computation. At that time, the control device 503 works as a control unit. A program executed with the general computation device 502 is stored in the main storage device 501 (storage unit). When the main storage device 501 has an insufficient memory capacity, the auxiliary storage device 504 that is similarly a storage unit is used. The input device 505 is used for inputting data, a program, and the like, and the output device 506 is used for outputting results. For the input device 505, not only a manual input device such as a keyboard but also an interface for a network connection may be used. In addition, the interface serves as the output device as well.

As described in the first to sixth embodiments, N spin variables s_j^z(t) and N effective magnetic field variables B_eff,j^z(t) are iteratively determined in the local-field response computation in the embodiments. The local-field response computing device 600 dedicatedly executes this iterative computation.

In the sixth embodiment, we performed the similar Processes 10a to 10d and obtained solutions s_j^zfd for the respective processes. The obtained solutions are transferred from the local-field response computing device 600 to the main storage device 501, and the energy H_p^q(t_m) and E_best are computed using the general computation device 502. That is, individual computations not belonging to the iterated computation is performed using the general computation device 502, and the dedication rate in the local-field response computing device 600 is increased.

Eighth Embodiment

High-difficulty problems in combinatorial optimization problems belong to NP-hard. In addition, problems classified into “P” and problems classified into “NP” both can be transferred to an NP-hard problem. Therefore, if an NP-hard combinatorial optimization problem is solved, almost all the combinatorial optimization problems are solved. A ground state-searching problems described Eq. (1) include NP-hard problems. In this embodiment, we show how to treat an NP-hard problem in accordance with Eq. (1) by using a maximum cut (MAX-CUT) problem that is a representative NP-hard problem as an example.

The maximum cut problem is a problem of graph theory. In the graph theory, a graph G is constructed from an node set V and a vertex set E, and is written as G=(V, E). An edge e is written as e={i, j} using two nodes. When a graph is defined by taking the direction of edges e into consideration, the graph is called a directed graph. When a graph is defined without taking the direction into consideration, the graph is called an undirected graph. For an edge e, a weight is also defined, and it is written as w_ij and w_ji. For an undirected graph, w_ij=w_ji. Let us think of dividing nodes in a weighted undirected graph G=(V, E) into two groups. The MAX-CUT problem is to find a division method to maximize a total sum of weights of cut edges in the division problem. Let G₁=(V₁, E₁) and G₂=(V₂, E₂) be the divided two undirected graphs. Then, the MAX-CUT problem is to maximize Eq. (12).

$\begin{array}{cc}w\left(V_1\right)=\sum _{i\in V_1,j\in V_2}w_{\mathrm{ij}}& \left(12\right)\end{array}$

If s_i=1 is set for edge i∈V₁, and s_j=−1 is set for edge j∈V₂, Eq. (13) is derived.

$\begin{array}{cc}w\left(V_1\right)=\sum _{i\in V_1,j\in V_2}w_{\mathrm{ij}}=\frac{1}{4}\sum _{\underset{\left(i\ne j\right)}{i,j\in V}}w_{\mathrm{ij}}s_is_j=\frac{1}{4}\sum _{\underset{\left(i\ne j\right)}{i,j\in V}}w_{\mathrm{ij}}-\frac{1}{2}\sum _{\underset{\left(i\ne j\right)}{i>j}}w_{\mathrm{ij}}s_is_j& \left(13\right)\end{array}$

Because the first term in the rightest side is a constant once a graph G is given, the MAX-CUT problem becomes a problem of minimizing Σ_i>jw_ijs_is_j. Because the Hamiltonian of the Ising spin glass model is given by Eq. (14),

$\begin{array}{cc}{\hat{H}}_p=-\sum _{i>j}J_{\mathrm{ij}}{\hat{\sigma }}_i^z{\hat{\sigma }}_j^z-\sum _jg_j{\hat{\sigma }}_j^z& \left(14\right)\end{array}$

the MAX-CUT problem is equivalent to the ground state-searching problem of Eq. (14) with J_ij=−w_ij and g_j=0.

Ninth Embodiment

In the second embodiment illustrated in FIG. 4, quantum entanglement is taken into consideration all over the variables at N sites corresponding to N spins in the computation. However, as the value of |s_j^z| increases, the probability with which a spin is inverted decreases. In this situation, spins may be assumed not to invert after a predetermined time, i.e., the effect of quantum entanglement may be ignored after the predetermined time. In other words, the computation at time t_i may be changed such that the condition of Step 113 in FIG. 4 is changed to l=x<N, and that if “Yes” at Step 113, the spins after x are not subjected to Step 112, and Steps 102a and 103 are performed individually and in parallel for all remaining sites. In this way, it is possible to shorten a process time while securing accuracy to some degree.

In this embodiment, the same function as that configured with software can be achieved with hardware such as an FPGA (Field Programmable Gate Array) and an ASIC (Application Specific Integrated Circuit).

The invention is not limited to the above embodiments, and various modifications can be made. For example, some configurations of a certain embodiment may be replaced with the configurations of another embodiment, and the configuration of the other embodiment may also be added to the configuration of a certain embodiment. In addition, part of the configurations of each embodiment may be added to or replaced with part of the configurations of other embodiments, and some of the configurations may be omitted.

Claims

1. A computing apparatus which includes a computing unit, a storage unit, and a control unit, and performs computation under the control of the control unit while transferring data between the storage unit and the computing unit, wherein N variables sjz (j=1, 2, . . . , N) take a range of −1≤sjz≤1, and an assignment is set with coefficients gj indicating local terms and coefficients Jkj (k, j=1, 2, . . . , N) indicating inter-variable interactions,time is divided into m, and the computing unit discretely performs computation from t=t0 (t0=0) to tm (tm≤τ),variables Beff,jz(ti) and sjz(ti) at each time ti (i=1, 2, . . . , m) are determined in this order, Beff,jz(ti) is a function of skz (ti−1) Jkj, gj, and ti, sjz(ti) is a function of Beff,jz (ti) and ti, and initial values at time t0 are set as Bjz(t0)=0 and sjz(t0)=0,for determining Beff,jz(ti) and sjz(ti) at time ti (i=1, 2, . . . , m), first, sjz(ti−1) are put in descending order such that |sm1z(ti−1)|≤sm2z(ti−1)|≤|sm3z(ti−1)|≤ . . . ≤|smNz(ti−1)|, then, Beff,m1z(ti) and sm1z(ti) at site m1 are determined at the first time, and sm1z(ti−1) is set to be sm1z(ti−1)=sgn(sm1z(ti))|sm1z(ti−1)|, next, Beff,m2z(ti) and sm2z(ti) at site m2 are determined and sm2z(ti−1) is set to be sm2z(ti−1)=sgn(sm2z(ti))|sm2z(ti−1)|, the computation at site m3 is performed similarly, and the computation up to site mx (herein, 3≤x≤N) is performed similarly for the computation at time ti, andthe variable sjz approaches −1 or 1 as the time step progresses from t=t0 to t=tm, and a solution is determined as sjzfd=−1 if sjz<0 and as sjzfd=1 if sjz>0.

2. The computing apparatus according to claim 1, wherein Beff,jz(ti) is determined by Beff,jz(ti)=(ti/τ)·(Σk(≠j)Jkjskz(ti−1)+gj).

3. The computing apparatus according to claim 1, wherein Beff,jx(ti)=γ(1−ti/τ) is set using a constant γ, θ is define by tan θ=Beff,jz(ti)/Beff,jx(ti), sjz(ti) is determined by sjz(ti)=sin θ, and thus sjz(ti)=f(Beff,jz(ti), ti)=sin{arctan(Beff,jz(ti)/Beff,jx(ti))} is obtained using a function f.

4. The computing apparatus according to claim 3, wherein correction parameters rs and rb are added to the function f,θ is defined by tan θ=rb·Beff,jz(ti)/Beff,jx(ti), and sjz(ti) is determined by sjz(ti)=rs·sin θ,and thus, the function f becomes f(Beff,jz(ti), ti)=rs·sin{arctan(rb·Beff,jz(ti)/Beff,jx(ti))}.

5. The computing apparatus according to claim 2, wherein ci=(Σk(skz(ti−1))2/N)1/2 and gjnorm(ti)=ci·gj are set, and Beff,jz(ti) is determined by Beff,jz(ti)=(ti/τ)·(Σk(≠j)Jkjsk(ti−1)+gjnorm(ti)).

6. The computing apparatus according to claim 5, wherein Beff,jz(ti) is determined by Beff,jz(ti)=(ti/τ)·(Σk(≠j)Jkjsk(ti−1)+ca·gjnorm(ti)) using a parameter ca.

7. The computing apparatus according to claim 4, wherein δrb≡1−rb is defined with respect to the correction parameter rb, and δrb is given as δrb(t)∝Σk(≠j)Jkj2.

8. The computing apparatus according to claim 5, wherein Bjz0(ti)=(Σk(≠j)Jkjskz(ti−1)+gjnorm(ti)) is defined, Bjz(ti)=(1−u)Bjz0(ti)+uBjz(ti−1) is defined using a parameter u satisfying 0≤u≤1, and Beff,jz(ti) is determined by Beff,jz(ti)=Bjz(ti)·ti/τ.

9. The computing apparatus according to claim 1, wherein the computation for determining sjzfd described in claim 1 is performed several times, a parameter div is set to a value as large as m, initial values at the second and subsequent computations are set as sjz(t0)=−sjzfd/div using the solution sjzfd for the last computation or are set as sjz(t0)=1/div or sjz(t0)=−1/div using a random number, Hp=−Σk>jJkjskzfd(ti)sjzfd−Σjgjsjzfd is calculated for each computation, and the final solution is sjzfd giving the minimum Hp in the repeated computations.

10. The computing apparatus according to claim 1, wherein after site mx, the computation of time ti is performed for all remaining sites independently and in parallel.

11. A computing method which uses a computing apparatus including a computing unit, a storage unit, and a control unit, and performs a computation under the control of the control unit while transferring data between the storage unit and the computing unit, wherein N variables sjz (j=1, 2, . . . , N) take a range of −1≤sjz≤1, and an assignment is set with coefficients gj indicating local terms and coefficients Jkj (k, j=1, 2, . . . , N) indicating inter-variable interactions,time is divided into m, and the computing unit discretely performs computation from t=t0 (t0=0) to tm (tm≤τ),variables Beff,jz(ti) and sjz(ti) at each time ti (i=1, 2, . . . , m) are determined in this order,Beff,jz(ti) is a function of skz(ti−1), Jkj, gj, and ti, sjz(ti) is a function of Beff,jz(ti) and ti, and initial values at time t0 are set as Bjz(t0)=0 and sjz(t0)=0,for determining Beff,jz(ti) and sjz(ti) at time ti (i=1, 2, . . . , m), first, sjz(ti−1) are put in descending order such that |sm1z(ti−1)|≤|sm2z(ti−1)|≤|sm3z(ti−1)|≤ . . . ≤|smNz(ti−1)|,Beff,m1z(ti) and sm1z(ti) at site m1 are determined at the first time, and sm1z(ti−1) is set to be sm1z(ti−1)=sgn(sm1z(ti))|sm1z(ti−1)|, for the following sites, the same computation is performed up to site mx (herein, 1≤x≤N) for the computation at time ti, andvariables sjz approach −1 or 1 as the time step progresses from t=t0 to t=tm, and a solution is determined as sjzfd=−1 if sjz<0 and as sjzfd=1 if sjz>0.

12. The computing method according to claim 11, wherein the same computation is performed up to site mN after site m1 for the computation at time ti.

13. The computing method according to claim 11, wherein after site mx, all remaining sites are processed independently and in parallel to perform the computation at time ti.