WO2016157333A1 - 計算機、及び演算プログラム - Google Patents

計算機、及び演算プログラム Download PDF

Info

Publication number
WO2016157333A1
WO2016157333A1 PCT/JP2015/059765 JP2015059765W WO2016157333A1 WO 2016157333 A1 WO2016157333 A1 WO 2016157333A1 JP 2015059765 W JP2015059765 W JP 2015059765W WO 2016157333 A1 WO2016157333 A1 WO 2016157333A1
Authority
WO
WIPO (PCT)
Prior art keywords
time
spin
variable
calculation
pina
Prior art date
Application number
PCT/JP2015/059765
Other languages
English (en)
French (fr)
Inventor
辰也 戸丸
Original Assignee
株式会社日立製作所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社日立製作所 filed Critical 株式会社日立製作所
Priority to US15/561,566 priority Critical patent/US10318607B2/en
Priority to PCT/JP2015/059765 priority patent/WO2016157333A1/ja
Priority to JP2017508862A priority patent/JP6291129B2/ja
Publication of WO2016157333A1 publication Critical patent/WO2016157333A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N10/00—Quantum computing, i.e. information processing based on quantum-mechanical phenomena

Definitions

  • the present invention relates to a computer that enables high-speed computation for inverse problems and combinatorial optimization problems that require exhaustive search.
  • Big data has many problems that require complex analysis. For example, when a certain result is obtained, there is a case where it is desired to find the cause. This is called an inverse problem. The more complicated the phenomenon is, the more difficult it is to find the cause, and there is generally no efficient algorithm for obtaining the initial value from the result. In the worst case, an exhaustive search must be performed for the initial value. This is one of the difficult problems in big data. Alternatively, there are many problems such as selecting an optimal solution from many options based on big data. In this case, if all possibilities are taken into account, an exhaustive search becomes necessary. against this background, there is a need for a computer that can efficiently solve problems that require exhaustive search.
  • a quantum computer consists of basic elements called qubits, and realizes "0" and "1" simultaneously. Therefore, all solution candidates can be calculated as initial values at the same time, and there is a possibility that an exhaustive search can be realized.
  • quantum computers need to maintain quantum coherence over the entire computation time, and there is no prospect of realizing this.
  • Non-Patent Document 1 a technique called adiabatic quantum computation has been attracting attention.
  • H ⁇ p be the Hamiltonian of the physical system that sets the problem.
  • the Hamiltonian is not H ⁇ p at the start of the calculation, but another Hamiltonian H ⁇ 0 whose ground state is clear and easy to prepare.
  • H ⁇ p the Hamiltonian of the physical system that sets the problem.
  • adiabatic quantum computation can be applied to problems that require exhaustive search, and reaches a solution in a one-way process.
  • the calculation process needs to follow the Schrödinger equation (2), it is necessary to maintain quantum coherence as in the case of the quantum computer.
  • quantum computers repeat the gate operation between 1 qubits or 2 qubits
  • adiabatic quantum computations interact all over the entire qubit system, and the concept of coherence is Different. For example, consider a gate operation to a certain qubit. At this time, if there is interaction between the qubit and other qubits, it causes decoherence, but in adiabatic quantum computation, all qubits interact simultaneously, so in this case Does not become decoherence. Reflecting this difference, adiabatic quantum computation is considered to be more robust against decoherence than quantum computers.
  • adiabatic quantum computation is effective for difficult problems that require exhaustive search.
  • it still requires quantum coherence, and when using a superconducting flux qubit, a cryogenic cooling device is also required. It is a problem to be solved to eliminate these two requirements and to provide a practical computer.
  • An object of the present invention is to solve the above-described problems and provide a computer and a calculation program that do not require quantum coherence or a cryogenic cooling device.
  • the spin is a variable in the operation, and the problem to be solved is set by the interaction between spins and the local field acting for each spin.
  • Each spin develops in time as the direction is determined according to the effective magnetic field determined by the external magnetic field at each site and all the interactions between the spins at time t.
  • the spin direction is not completely aligned with the effective magnetic field, but the direction is corrected in a quantum mechanical manner so that the system substantially maintains the ground state.
  • the vibration related to the direction of the spin in the time evolution is suppressed, and the convergence of the solution is improved.
  • FIG. 20 is a block diagram illustrating a configuration example of a computer according to an eighth embodiment.
  • FIG. 20 is a block diagram illustrating another configuration example of the computer according to the eighth embodiment. It is a block diagram which shows the example of 1 structure of the local field response calculating apparatus contained in the computer based on Example 9.
  • FIG. It is a block diagram which shows one structural example of the optical part of FIG.
  • FIG. It is a block diagram which shows the example of 1 structure of the local field response calculating apparatus contained in the computer based on Example 10.
  • FIG. It is a block diagram which shows the example of 1 structure of the local field response calculating apparatus contained in the computer based on Example 11.
  • FIG. It is a block diagram which shows the example of 1 structure of the local field response calculating apparatus contained in the computer based on Example 12.
  • FIG. 12 shows the example of 1 structure of the local field response calculating apparatus contained in the computer based on Example 12.
  • notations such as “first”, “second”, and “third” are attached to identify the constituent elements, and do not necessarily limit the number or order.
  • a number for identifying a component is used for each context, and a number used in one context does not necessarily indicate the same configuration in another context. Further, it does not preclude that a component identified by a certain number also functions as a component identified by another number.
  • Adiabatic quantum computation also known as quantum annealing, is an evolution of the classic annealing concept to quantum mechanics.
  • the adiabatic quantum computation is inherently capable of classical operation, and can be interpreted as a quantum mechanical effect added to improve performance in terms of high speed and correct solution rate. Therefore, in the present invention, the arithmetic unit itself is classical, and by introducing a parameter determined in a quantum mechanical manner into the arithmetic process, a classical arithmetic method and apparatus including a quantum mechanical effect is realized.
  • the following embodiment describes a classical algorithm for obtaining a ground state as a solution while explaining the relevance to adiabatic quantum computation, and a device for realizing it.
  • a typical form described in the following embodiments is a computer that includes a calculation unit, a storage unit, and a control unit, and performs calculations while exchanging data between the storage unit and the calculation unit under the control of the control unit.
  • the coefficient g pinb is, for example, a value from 50% to 200% of the average value of
  • the correction term ⁇ g j ′ is added to g j ′ only for a certain site j ′, and the size of g j ′ is increased only for the site j ′. it can.
  • the correction term ⁇ g j ′ is, for example, a value of 10% to 100% of the average value of
  • Example 1 the principle of the present invention is described through a transition from the quantum mechanical description to the classical form.
  • J ij and g j are task setting parameters, and ⁇ ⁇ j z is the z component of Pauli's spin matrix and takes an eigenvalue of ⁇ 1.
  • i and j represent spin sites.
  • Ising spin is a variable that can take only ⁇ 1 as a value.
  • the eigenvalue of ⁇ ⁇ j z is ⁇ 1, so it is an Ising spin system.
  • the Ising spin in equation (3) does not have to be a literal spin, and can be physically anything as long as the Hamiltonian is described by equation (3).
  • the logic circuit high and low can be associated with ⁇ 1
  • the longitudinal and transverse polarization of light can be associated with ⁇ 1
  • the phase of 0, ⁇ can be associated with ⁇ 1. is there.
  • Equation (4) corresponds to the application of a transverse magnetic field, and the ground state is when all spins are in the x direction ( ⁇ > 0).
  • the problem setting Hamiltonian was defined as an Ising spin system with only the z component, but the x component of the spin appears in Equation (4). Therefore, the spin in the calculation process is not Ising but vector-like (Bloch vector).
  • t 0, the start is made with the Hamiltonian in the equation (4), but the Hamiltonian is gradually changed with the progress of the time t, and finally the Hamiltonian described in the equation (3) is used to obtain the ground state as a solution.
  • ⁇ ⁇ represents the three components of Pauli's spin matrix as a vector.
  • the spin direction can be defined by ⁇ ⁇ j z > / ⁇ ⁇ j x >, if the spin direction follows the effective magnetic field, the spin direction is determined by equation (7).
  • Equation (7) is a quantum mechanical description but has an expected value, so it is a relational equation for classical quantities, unlike equations (1)-(6).
  • equations (1)-(6) there is no non-local correlation (quantum drift) of quantum mechanics, so the spin direction should be completely determined by the local field at each site, and Equation (7) determines the behavior of the classical spin system. Since there is a nonlocal correlation in the quantum system, Equation (7) is modified, but this will be described in Example 2 and below, and in this example, Equation (7) is used to describe the basic form of the invention. Describes the classical system determined by).
  • FIG. 1 shows a timing chart (procedure 100) for obtaining the ground state of the spin system.
  • the spin of site j is represented by s j instead of ⁇ ⁇ j .
  • the effective magnetic field B j in FIG. 1 is a classical quantity.
  • a right effective magnetic field B j is applied at all sites, and all spins s j are initialized to the right.
  • the magnetic field in the z-axis direction and the spin-to-spin interaction are gradually applied.
  • the time t is continuous, but it can be made discrete in the actual calculation process to improve convenience. The discrete case is described below.
  • the spin according to the present invention is a vector spin because not only the z component but also the x component is added.
  • the behavior as a vector can also be understood from FIG.
  • a three-dimensional vector having a size of 1 this is called a Bloch vector and a state can be described by a point on the sphere
  • only two dimensions are used. (The state can be described by a point on the circle).
  • the two-dimensional spin vector can be described only by a semicircle. If s j z is specified by [ ⁇ 1, 1], the two-dimensional spin vector is determined by one variable of s j z . Therefore, although the spin of the present invention is a two-dimensional vector, it can also be expressed as a one-dimensional continuous variable with a range of [ ⁇ 1, 1].
  • Equation (8) is a rewrite of Equation (7) into a notation related to classical quantities, it does not have a ⁇ •> symbol.
  • the magnitude of the spin vector is 1.
  • FIG. 2 shows a summary of the above algorithm in a flowchart.
  • t m ⁇ .
  • variable interaction J ij The strength of the correlation between the temperature and the activities of public facilities and private houses is expressed through the variable interaction J ij .
  • ⁇ ⁇ j z is a variable representing power distribution, but correlates with the movement of people and the opening of public facilities. Therefore, depending on the solution obtained, it can be interpreted as “a certain public facility should be closed”.
  • the above is a simple example of expressing a specific problem with Expression (3).
  • the specific problems to which this embodiment can be applied are not limited to the power supply management problems exemplified above, but include travel route optimization, vehicle guidance for avoiding traffic congestion, circuit design, product supply management, scheduling, financial asset selection, etc. It can be applied to solve many problems.
  • Example 1 the expected value was taken based on the quantum mechanical equation to shift to the classical quantity, and the algorithm based on the classical quantity was explained using FIG. 1 and FIG. Since Example 1 was mainly intended to explain the algorithm of the present invention, it was described without including quantum mechanical effects. However, if the quantum mechanical effect is added, it can be expected to improve the accuracy rate and the calculation speed. Therefore, in the second embodiment, although the algorithm itself is classic, a method of adding a correction parameter based on quantum mechanics to improve performance will be described.
  • the features of quantum mechanics include linear superposition and quantum distortion (nonlocal correlation). For example, consider a qubit that can take two states,
  • the linear superposition state is a sum of two states such as
  • ⁇ > ⁇
  • 0>, and if s j z (t k ) ⁇ 1, the state is
  • 0> i
  • 1> ⁇ i
  • ⁇ > i ⁇
  • the magnitude of the spin vector is not preserved to 1 when there is quantum distortion although it is an example.
  • the magnitude of the spin vector should be a constant 1, but if there is a quantum droop, the spin vector will not have a magnitude of 1.
  • s j z (t with ⁇ defined as tan ⁇ ⁇ B j z (t)> / ⁇ B j x (t)> as a parameter.
  • this method does not reflect the nature of the quantum drowning inherent in this system. Therefore, let us consider reflecting the nature of quantum drowning.
  • the correction parameters r s and r B originate from quantum distortion and are preferably finely controlled depending on t k , s j z (t k ), and s j x (t k ).
  • r B is an amount that can change the sign and reflects the quantum blur best.
  • r s and r B are not given site dependency, but if the characteristics of each site are known in advance, r s and r B are set accordingly. If you make it site-dependent, you can expect an improvement in the correct answer rate.
  • FIG. 3 shows a flowchart when r s and r B are introduced. The difference from the flowchart of FIG. 2 is that steps 103, 105, and 107 are changed to steps 103a, 105a, and 107a including correction parameters r s and r B , respectively.
  • Equation (12A) the relaxation term (pinning term) is added to Equation (10) to make the effective magnetic field into Equation (12A).
  • the third term is a relaxation term.
  • the relaxation term is an additional term for improving the convergence of the solution and needs to be sufficiently smaller than
  • Relaxation term (third term) of the formula (12A) does not depend on the size of depends only on the sign of s j z s j z.
  • the first term depends on the size of s i z . Therefore, there is also a method of making the third term depend on the size of s j z . In that case, Expression (12B) is obtained.
  • the coefficient g pinb is typically about the average value of
  • s j z ⁇ 0.
  • the relaxation term is also useful for this induction.
  • one site referred to as j site
  • s j z (t 0 ) 1
  • the orientation of other sites is determined at the initial stage based on the sign of J ij with reference to j site. To do. In this way, an appropriate zeroth approximation is set, and then converges to one of the correct answers through subsequent local field response calculations.
  • the relaxation term contributes as follows.
  • FIG. 4 shows a flowchart when a relaxation term is added.
  • 4A corresponds to the equation (12A)
  • FIG. 4B corresponds to the equation (12B).
  • steps 102, 104, and 106 are changed to steps 102b, 104b, and 106b including a relaxation term.
  • ⁇ g j ′ is added to the local field term g j ′ in equations (12A) and (12B) to obtain g j ′ ⁇ g j ′ + ⁇ g j ′ .
  • the relaxation term improves the convergence of the solution and converges to one of the solutions even when the number of degeneracy is large.
  • the relaxation term improves the convergence of the solution and converges to one of the solutions even when the number of degeneracy is large.
  • formulas (12A) and (12B) only the local field term ⁇ g j ′ is added to only one site j ′ to strongly induce only the spin at the j ′ site in a certain direction. This leads to one of the degenerate solutions based on the j 'site.
  • the convergence of the solution is improved by introducing a relaxation term.
  • FIG. 5 shows a flowchart of an example of the algorithm in this case. 5A-5D correspond to FIGS. 4A-4D, respectively.
  • FIG. 6 is a flowchart illustrating an example of an algorithm related to a final solution determination method.
  • ⁇ ⁇ j z in Equation (3) has an eigenvalue of ⁇ 1, it is determined at each stage of the calculation process whether the eigenvalue of ⁇ ⁇ j z is 1 or -1 according to the sign of s j z .
  • the fourth method terminates the operation when all s j z converge as in the second method.
  • the final solution is not determined from s j z at the time of termination, but the energy at each time is calculated in the same way as the third method, and s j z at the time when the lowest energy is given ( The final solution is determined from t k ′ ) (125). The user decides which method to use.
  • the second method evaluates the possibility of spin reversal at each time and sets a time interval based on the evaluation. For example: If the magnitude of
  • a specific example of the time interval determination method is as follows. Let the minimum time interval be ⁇ t min .
  • the minimum value of the time interval is ⁇ t min and the maximum value is 100 ⁇ ⁇ t min .
  • the user decides which method to use.
  • Example 1-7 the calculation principle and calculation algorithm have been described.
  • a configuration example of a computer that operates this algorithm as a program will be described first.
  • FIG. 7 shows an example of the computer configuration of this embodiment.
  • FIG. 7 has basically the same configuration as that of an ordinary computer. Data is transferred from the main storage device 201 serving as a storage unit to the arithmetic device 202 serving as a calculation unit, and the data is returned to the main storage device 201 after the calculation. The operation is repeated by repeating this.
  • the control tower at that time is a control device 203 as a control unit.
  • the calculation in the configuration of the present embodiment is executed by the calculation device 202 for each time and for each site according to the flow of FIGS.
  • the program executed by the arithmetic device 202 is stored in the main storage device 201 which is a storage unit. If the main storage device 201 has insufficient storage capacity, the auxiliary storage device 204, which is also a storage unit, is used.
  • An input device 205 is used to input data and programs, and an output device 206 is used to output the results.
  • the input device 205 includes an interface for network connection in addition to a manual input device such as a keyboard. This interface also serves as an output device. 7 applies the algorithm described in Embodiment 1-7 as a program, but the apparatus itself uses a normal computer.
  • Fig. 8 shows a configuration diagram that achieves this.
  • a local field response calculation device 1000 is included as a calculation unit. This is a difference in comparison with the configuration of FIG.
  • the local field response calculation device 1000 is a dedicated calculation unit for executing the above-described algorithm, and performs only the local field response calculation described in the embodiment 1-7, and general processing is performed by the calculation device 202.
  • Information used in the local field response calculation device 1000 is transferred from the main storage device 201.
  • Necessary information includes task setting parameters such as interaction J ij between variables and local field g j , time parameter t k , correction parameters r s and r B related to quantum distortion , and coefficients g pina and g pinb of relaxation terms Etc. Processing such as synchronization is the same as that of the computer having the configuration shown in FIG.
  • the transfer value from the local field response arithmetic unit 1000 in the intermediate stage of the calculation is used for, for example, calculating the Ising spin Hamiltonian energy in the intermediate stage of the calculation and evaluating s min (t k ) 2 and s ave (t k ) 2 . Use.
  • the final solution is determined using the energy value of the Ising spin Hamiltonian at the intermediate stage. If it is determined whether s j zd is +1 or ⁇ 1, the calculation of the energy value of the Ising spin Hamiltonian is simple, and a normal arithmetic unit 202 is used.
  • the local field response calculation device 1000 is a dedicated calculation unit for the local field response calculation described in the embodiment 1-4, and a normal calculation device 202 is used for other processing.
  • the local field response calculation apparatus 1000 described in the eighth embodiment can be realized by various methods. In this embodiment, first, a method for effectively using the parallelism of light will be described.
  • FIG. 9 shows an overall view.
  • Information such as ⁇ , g j , and g pina required for calculation is sent from the main storage device 201 to the control unit 1100, and information about J ij is sent to the variable mask 1120.
  • the output intensity of the LED array 1110 reflects the value of the variable s j z .
  • phase information is not used. Therefore, an incoherent light source such as an LED array is used as the light source 1110.
  • the measurement time in the detector array 1130 is sufficiently long so that the interference effect between the output lights of the LD array does not appear.
  • the output light from the LED array 1110 spreads only in the horizontal direction, attenuates the amount of light according to J ij in the variable mask 1120, and this time converges only in the vertical direction and converts it into an electrical signal by the detector array 1130.
  • FIG. 10 is a block diagram showing a configuration example of the optical unit in FIG.
  • the optical system from the LED array 1110 to the detector array 1130 uses, for example, a lens system as shown in FIG.
  • , and the difference between the two outputs s j z s j z + ⁇
  • s j z + + s
  • Two detector arrays 1130 are also paired in conjunction with the light source side. As a result, it is possible to cope with the variable mask 1120 that cannot take a negative value.
  • the signal to be obtained by the detector array 1130 is b j z ⁇ i J ij s i z .
  • Each of the detector array 1130 and two pair b j z + ⁇ i ( J ij + s j z + + J ij - s j z-) and
  • ⁇ i (J ij +
  • s j z + ), and the difference between the detectors is detected b j z b j z + ⁇
  • b j z + + b j z -Becomes a signal.
  • B j z can be obtained based on equations (10), (12A), and (12B) by adding the terms g j and g pina .
  • This calculation is performed by the control unit 1100.
  • the calculation to obtain the s j z from B j z also performed by the control unit 1100, sends the value of s j z to the LED array 1110.
  • the control unit 1100 repeats the same processing, and is a dedicated circuit for that purpose. Further, s j z at each time is transferred to the main storage device 201 and used for analysis.
  • the optical system from the LED array 1110 to the detector array 1130 was free space. This part can also be realized by an optical system using a waveguide.
  • FIG. 11 shows such a case.
  • the output light from the LED array (LD array) 1110 is divided by a demultiplexer 1115, passes through a variable attenuator 1120 (variable mask in FIG. 9) whose transmittance is set based on J ij, and is collected by a multiplexer 1125.
  • the light is received and received by the detector array 1130. Only the spatial optical system in FIG. 9 is changed to a waveguide optical system, and the operation principle is the same. Accordingly, the LED array (LD array) 1110 represents s j z in pairs, and the detector array 1130 also operates in pairs.
  • Example 9 s j z was expressed by combining two light sources. If s j z + and s j z ⁇ are expressed using polarized light, one light source can be used for s j z .
  • FIG. 12 shows such a case.
  • Each LD in the light source 1111 is for representing s j z .
  • Acos ( ⁇ / 4 + ⁇ ).
  • the polarization-modulated light is separated by the polarization separator 1113 and guided to the branching filter 1115.
  • the operation from the duplexer 1115 to the detector 1130 is the same as that in FIG. Since the detector 1130 operates in pairs as in FIG.
  • the local field response calculation device 1000 can be realized not only by a method using light as in the ninth to eleventh embodiments but also by an electric circuit.
  • FIG. 13 shows a configuration example in that case.
  • Information on each spin s i z is temporarily held in the buffer array 1210.
  • Information on J ij is stored in the memory 1221.
  • g pinb is held at 1221
  • J ii g pinb .
  • the cells for each s i z in the buffer arrays 1210 and 1230 are composed of multiple bits, and are assumed to be a pseudo continuous quantity.
  • the influence of temperature in the present invention is estimated as follows. Bit manipulation is performed by the LED (LD) array 1110, the polarization modulator 1112, and the buffer arrays 1210 and 1230.
  • the present invention is not limited to the above-described embodiment, and includes various modifications.
  • a part of the configuration of one embodiment can be replaced with the configuration of another embodiment, and the configuration of another embodiment can be added to the configuration of one embodiment.

Landscapes

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

Abstract

全数探索を必要とするような課題に対して、量子コヒーレンスや極低温冷却装置を必要としない計算機、及び演算プログラムを提供する。　変数としてのスピンsjを局所的な有効磁場Bjに追従させることにより、問題設定系の基底状態に系を導く。その基底状態が解である。t ＝ 0において全サイトでx軸方向に磁場Bjが印加され、全スピンsjがx軸方向に初期化される。時間tの経過に従い、徐々にz軸方向の磁場とスピン間相互作用が加えられ、最終的にスピンは＋z方向あるいは－z方向となって、スピンsjのz成分がsj z ＝ ＋1あるいは－1となる。また、スピンsjの向きを有効磁場Bj方向に追従させる際に、スピンsjの向きを維持しようとする緩和項を導入して解の収束性を向上させる。

Description

計算機、及び演算プログラム
　本発明は、全数探索を必要とするような逆問題や組み合わせ最適化問題に対して高速演算を可能にする計算機に関するものである。
　ビッグデータといった言葉に代表されるように現代はデータが溢れている。情報科学では、この巨大データを如何に解析し如何に扱うかが最重要課題の一つになっている。ビッグデータでは複雑な解析を必要とする問題が少なくない。例えば、ある結果が得られた際にその原因を探りたい場合がある。これを逆問題と呼ぶことにする。現象が複雑であればあるほど原因を探るのは困難であり、一般に結果から初期値を求める効率的なアルゴリズムは存在しない。最悪の場合は初期値に関して全数探索しなければならない。これがビッグデータにおける困難問題の一つである。あるいはまた、ビッグデータを元に多くの選択肢の中から最適解を選ぶといった問題も多く存在する。この場合もすべての可能性を考慮するならば全数探索の必要性が出てくる。こういった背景から、全数探索を必要とする問題を効率的に解けるコンピュータが求められることになる。
　全数探索問題に対しては量子コンピュータへの期待が大きい。量子コンピュータは量子ビットと呼ばれる基本素子からなり“0”と“1”を同時に実現する。そのためすべての解候補を初期値として同時に計算可能で、まさに全数探索を実現しうる可能性を持っている。しかし、量子コンピュータは全計算時間に亘って量子コヒーレンスを維持する必要があり、これを実現する目処は立っていない。
　こういった中で注目されるようになってきたのが断熱量子計算と呼ばれる手法である（非特許文献１）。この方法は、ある物理系の基底状態が解になるように問題を変換し、基底状態を見つけることを通して解を得ようとするものである。問題を設定した物理系のハミルトニアンをH^pとする。但し、演算開始時点ではハミルトニアンをH^pとするのではなく、それとは別に基底状態が明確で準備しやすい別のハミルトニアンH^0とする。次に十分に時間を掛けてハミルトニアンをH^0からH^pに移行させる。十分に時間を掛ければ系は基底状態に居続け、ハミルトニアンH^pの基底状態が得られる。これが断熱量子計算の原理である。計算時間をτとすればハミルトニアンは式(1)となり
Figure JPOXMLDOC01-appb-M000001
　式(2)のシュレディンガー方程式に基づいて時間発展させて解を得る。
Figure JPOXMLDOC01-appb-M000002
　断熱量子計算は全数探索を必要とする問題に対しても適用可能で、一方向性の過程で解に到達する。しかし、計算過程が式(2)のシュレディンガー方程式に従う必要があるならば、量子コンピュータと同様に量子コヒーレンスの維持が必要になる。但し、量子コンピュータが1 量子ビットあるいは2 量子ビット間に対するゲート操作を繰り返すものであるのに対して、断熱量子計算は量子ビット系全体に亘って一斉に相互作用させるものであり、コヒーレンスの考え方が異なる。例えば、ある量子ビットへのゲート動作を考えてみる。この時もしその量子ビットと他の量子ビットとで相互作用があれば、それはディコヒーレンスの原因になるが、断熱量子計算ではすべての量子ビットを同時に相互作用させるので、この例のような場合にはディコヒーレンスにならない。この違いを反映して断熱量子計算は量子コンピュータに比べてディコヒーレンスに対して頑強であると考えられている。
　しかし、断熱量子計算にも課題がある。ディコヒーレンスに関して量子コンピュータと比べて頑強になったとしても演算過程が式(2)のシュレディンガー方程式に従うならばやはり十分なコヒーレンスが必要である。また断熱量子計算を実現しうる系が超伝導磁束量子ビット系であるのも課題になる（特許文献１、非特許文献２）。なぜなら、超伝導を用いる場合には極低温冷却装置を必要とするからである。極低温を必要とすることは実用的なコンピュータ実現のためには課題になる。
特表２００９－５２４８５７号公報
E. Farhi, et al., "A quantum adiabatic evolution algorithm applied to random instances of an NP－complete problem," Science 292, 472 (2001). A. P.－Ortiz, "Finding low－energy conformations of lattice protein models by quantum annealing," Scientific Reports 2, 571 (2012). F. Barahona, "On the computational complexity of Ising spin glass models," J. Phys. A: Math. Gen. 15, 3241 (1982).
　以上述べたように、断熱量子計算は全数探索を必要とするような難問に対して有効である。しかし、依然として量子コヒーレンスを必要とし、また超伝導磁束量子ビットを用いる場合には極低温冷却装置も必要である。これら２つの必要条件を無くし、実用的な計算機を提供することが解決すべき課題である。
　本発明の目的は、上記の課題を解決し、量子コヒーレンスや極低温冷却装置を必要としない計算機、及び演算プログラムを提供することにある。
　スピンを演算における変数とし、解こうとする問題をスピン間相互作用とスピンごとに作用する局所場で設定する。時刻t＝0において外部磁場により全スピンを一方向に向かせ、時刻t＝τで外部磁場がゼロになるように外部磁場を徐々に小さくする。各スピンは時刻tにおける各サイトの外部磁場及びスピン間相互作用のすべての作用で決まる有効磁場に従い向きが定まるとして時間発展させる。その際、スピンの向きが有効磁場に完全に揃うのではなく、量子力学的に補正された向きとすることにより、系が基底状態をほぼ維持するようにする。
　また、時間発展の際に各スピンを元の向きに維持する項（緩和項）を有効磁場に加え、解の収束性を向上させる。
　本方法は量子力学的補正を加えるものの古典系で動作するものである。そのため、量子コヒーレンスを考慮する必要がなくなり使用可能なリソースが広範囲になる。ビットに関するエネルギースケールを温度のエネルギースケールよりも十分に大きくすれば、温度揺らぎも無視可能で、極低温装置のような特殊装置・特殊環境が不要になる。
　また、緩和項を加えたことにより、時間発展におけるスピンの向きに関する振動が抑えられ、解の収束性が向上する。
本発明のアルゴリズムを原理的に説明する模式図である。 実施例１に係る、アルゴリズムの一例をフローチャートとして示す図である。 実施例２に係る、アルゴリズムの一例をフローチャートとして示す図である。 実施例３に係る、アルゴリズムの一例をフローチャートとして示す図である。 実施例３に係る、アルゴリズムの他の例をフローチャートとして示す図である。 実施例３に係る、アルゴリズムの他の例をフローチャートとして示す図である。 実施例３に係る、アルゴリズムの他の例をフローチャートとして示す図である。 実施例４に係る、アルゴリズムの一例をフローチャートとして示す図である。 実施例４に係る、アルゴリズムの他の例をフローチャートとして示す図である。 実施例４に係る、アルゴリズムの他の例をフローチャートとして示す図である。 実施例４に係る、アルゴリズムの他の例をフローチャートとして示す図である。 実施例５に係る、最終的な解決定法に関するアルゴリズムの一例をフローチャートとして示す図である。 実施例５に係る、最終的な解決定法に関するアルゴリズムの他の例をフローチャートとして示す図である。 実施例５に係る、最終的な解決定法に関するアルゴリズムの他の例をフローチャートとして示す図である。 実施例５に係る、最終的な解決定法に関するアルゴリズムの他の例をフローチャートとして示す図である。 実施例８に係る、計算機の一構成例を示すブロック図である。 実施例８に係る、計算機の他の構成例を示すブロック図である。 実施例９に係る、計算機に含まれる局所場応答演算装置の一構成例を示す構成図である。 図９の光学部の一構成例を示す構成図である。 実施例１０に係る、計算機に含まれる局所場応答演算装置の一構成例を示すブロック図である。 実施例１１に係る、計算機に含まれる局所場応答演算装置の一構成例を示すブロック図である。 実施例１２に係る、計算機に含まれる局所場応答演算装置の一構成例を示すブロック図である。
　以下、演算の原理も交えて、図面に従い本発明の各種実施例を説明する。ただし、本発明は以下に示す実施の形態の記載内容に限定して解釈されるものではない。本発明の思想ないし趣旨から逸脱しない範囲で、その具体的構成を変更し得ることは当業者であれば容易に理解される。
　以下に説明する発明の構成において、同一部分又は同様な機能を有する部分には同一の符号を異なる図面間で共通して用い、重複する説明は省略することがある。
　本明細書等における「第１」、「第２」、「第３」などの表記は、構成要素を識別するために付するものであり、必ずしも、数または順序を限定するものではない。また、構成要素の識別のための番号は文脈毎に用いられ、一つの文脈で用いた番号が、他の文脈で必ずしも同一の構成を示すとは限らない。また、ある番号で識別された構成要素が、他の番号で識別された構成要素の機能を兼ねることを妨げるものではない。
　図面等において示す各構成の位置、大きさ、形状、範囲などは、発明の理解を容易にするため、実際の位置、大きさ、形状、範囲などを表していない場合がある。このため、本発明は、必ずしも、図面等に開示された位置、大きさ、形状、範囲などに限定されない。
　断熱量子計算は別名で量子アニールとも呼ばれ、古典的な焼きなましの概念を量子力学に発展させたものである。即ち、断熱量子計算は本来古典的動作が可能で、高速性や解の正解率に関して性能を向上させるために量子力学的効果が付加されたものとも解釈できる。そこで本発明では、演算器そのものは古典的とし、演算過程に量子力学的に定まるパラメタを導入することにより、古典的であるが量子力学的な効果を含んだ演算方法・装置を実現する。
　以上の概念に基づき、以下の実施例では断熱量子計算との関連性を説明しながら解としての基底状態を得る古典的アルゴリズムと、それを実現するための装置に関して述べる。
　以下の実施例で述べられる典型的な形態は、演算部、記憶部、制御部を具備し、制御部の制御により、記憶部と演算部との間でデータをやり取りしながら演算を行う計算機であって、N個の変数sj z(j ＝ 1, 2, …, N)が-1 ≦ sj z ≦ 1の値域を取り、局所場gjと変数間相互作用Jij (i, j ＝ 1, 2, …, N)によって課題の設定がなされる。演算部では、時刻をm分割して離散的にt ＝ t0 (t0 ＝ 0)からtm (tm ＝ τ)まで演算するものとし、各時刻tkにおける変数sj z(tk)を求めるに当たり、前時刻tk-1の変数si z(tk-1) (i ＝ 1, 2, .., N)の値と緩和項の係数gpinaあるいはgpinbを用いてBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ sgn(sj z(tk-1))・gpina}・tk/τあるいはBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ gpinb・sj z(tk-1)}・tk/τを求め、 前記変数sj z(tk)の値域が-1 ≦ sj z(tk) ≦ 1になるように関数fを定めてsj z(tk) ＝ f(Bj z(tk) , tk)とし、時刻ステップをt ＝ t0からt ＝ tmに進めるにつれて前記変数sj zを-1あるいは1に近づけ、最終的にsj z ＜ 0ならばsj zd ＝ -1、sj z ＞ 0ならばsj zd ＝ 1として解を定める。
　係数gpinbは、例えば|Jij|の平均値の50%から200%の値である。また、課題設定の局所場gjに関して、あるサイトj’に対してのみ補正項δgj’をgj’に加え、該サイトj’に対してのみgj’の大きさを大きくすることもできる。また、補正項δgj’は、例えば|Jij|の平均値の10%から100%の値である。
　実施例１では、量子力学的な記述から出発して古典的な形式に移行することを通して本発明の原理を述べる。
　式(3)で与えられるイジングスピン・ハミルトニアンの基底状態探索問題はNP困難と呼ばれる分類の問題を含み、有用な問題であることが知られている（非特許文献３）。
Figure JPOXMLDOC01-appb-M000003
　Jij及びgjが課題設定パラメタであり、σ^j zはパウリのスピン行列のz成分で±1の固有値を取る。i, jはスピンのサイトを表す。イジングスピンとは値として±1だけを取りうる変数のことで、式(3)ではσ^j zの固有値が±1であることによりイジングスピン系となっている。式(3)のイジングスピンは文字通りのスピンである必要はなく、ハミルトニアンが式(3)で記述されるのであれば物理的には何でも良い。例えば、ロジック回路のhighとlowを±1に対応付けることも可能であるし、光の縦偏波と横偏波を±1に対応付けることや0, πの位相を±1に対応付けることも可能である。本実施例の方法では断熱量子計算と同様に、時刻t ＝ 0において式(4)で与えられるハミルトニアンの基底状態に演算系を準備する。
Figure JPOXMLDOC01-appb-M000004
　γは全サイトjに一様に掛かる外場の大きさで決まる比例定数であり、σ^j xはパウリのスピン行列のx成分である。演算系がスピンそのものであれば外場とは磁場を意味する。式(4)は横磁場を印加したことに相当し、すべてのスピンがx方向を向いた場合(γ>0)が基底状態である。問題設定のハミルトニアンはz成分のみのイジングスピン系として定義されたが、式(4)にはスピンのx成分が登場している。従って、演算過程でのスピンはイジングではなくベクトル的（ブロッホベクトル）である。t ＝ 0では式(4)のハミルトニアンでスタートしたが、時刻tの進行と共に徐々にハミルトニアンを変化させ、最終的には式(3)で記述されるハミルトニアンにしてその基底状態を解として得る。
　スピンが外場に対してどのように応答するかを、１スピン系の場合にまず考察する。１スピン系のハミルトニアンは式(5)で与えられる。
Figure JPOXMLDOC01-appb-M000005
　ここでσ^はパウリのスピン行列の３成分をベクトルとして表示している。基底状態はスピンが磁場方向を向いた場合で、<・>を量子力学的期待値として<σ^> ＝ B/|B|と書ける。断熱過程では常に基底状態を維持しようとするのでスピンの向きは常に磁場の向きに追従する。
　以上の議論は多スピン系にも拡張できる。t ＝ 0ではハミルトニアンが式(4)で与えられる。これは全スピンに対して一様に磁場Bj x ＝γが印加されたことを意味する。t > 0では磁場のx成分が徐々に弱まりBj x ＝ γ(1 － t/τ)である。z成分に関してはスピン間相互作用があるために有効磁場としては式(6)になる。
Figure JPOXMLDOC01-appb-M000006
　スピンの向きは<σ^j z>/<σ^j x>で規定できるので、スピンの向きが有効磁場に追従するならば式(7)によりスピンの向きが定まる。
Figure JPOXMLDOC01-appb-M000007
　式(7)は量子力学的記述であるが期待値を取っているので、式(1) － (6)とは異なり古典量に関する関係式である。古典系では量子力学の非局所相関（量子縺れ）がないので、スピンの向きはサイトごとの局所場により完全に決まるはずであり、式(7)が古典的スピン系の振る舞いを決定する。量子系では非局所相関があるために式(7)は変形されることになるが、それに関しては実施例２以下で述べることとし、本実施例では発明の基本形態を述べるために式(7)で定まる古典系について記述する。
　図１にスピン系の基底状態を得るためのタイミングチャート（手順１００）を示す。図１の記述は古典量に関するものなので、サイトjのスピンをσ^jではなくsjにより表した。またそれに伴い、図１の有効磁場Bjは古典量である。t ＝ 0において全サイトで右向きの有効磁場Bjが印加され、全スピンsjが右向きに初期化される。時間tの経過に従い、徐々にz軸方向の磁場とスピン間相互作用が加えられ、最終的にスピンは＋z方向あるいは－z方向となって、スピンsjのz成分がsj z ＝ ＋1あるいは－1となる。時間tは連続的であることが理想であるが、実際の演算過程では離散的にして利便性を向上させることもできる。以下では離散的な場合を述べる。
　本発明に係るスピンはz成分だけでなくx成分が加わっているためにベクトル的なスピンになっている。図１からもベクトルとしての振る舞いが理解できる。ここまでy成分が登場してこなかったが、それは外場方向をxz面に取ったために外場のy成分が存在せず、従って<σ^y> ＝ 0となるためである。演算系のスピンとしては大きさ１の３次元ベクトル（これをブロッホベクトルと呼び、球面上の点で状態を記述できる）を想定しているが、本実施例の軸の取り方では２次元のみを考慮すればよい（円上の点で状態を記述できる）。またγは一定なのでBj x(t) > 0 (γ> 0)あるいはBj x(t) ＜ 0 (γ＜ 0)が成り立つ。この場合、２次元スピンベクトルは半円のみで記述できることになり、[－1, 1]でsj zを指定すればsj zの１変数で２次元スピンベクトルが定まる。従って、本発明のスピンは２次元ベクトルであるが、値域を[－1, 1]とする１次元連続変数として表記することもできる。
　図１の手順１００では時刻t ＝ tkにおいてサイトごとに有効磁場を求め、その値を用いて式(8)によりt ＝ tkにおけるスピンの向きを求める。
Figure JPOXMLDOC01-appb-M000008
　式(8)は式(7)を古典量に関する表記に書き改めたものなので<・>の記号が付いていない。
　次に、t ＝ tk＋1の有効磁場をt ＝ tkにおけるスピンの値を用いて求める。各時刻の有効磁場を具体的に書けば式(9)及び(10)となる。
Figure JPOXMLDOC01-appb-M000009
Figure JPOXMLDOC01-appb-M000010

　以下、図１の手順100で模式的に示した手順に従い、スピンと有効磁場を交互に求めていく。
　古典系ではスピンベクトルの大きさは１である。この場合スピンベクトルの各成分は、tanθ ＝ Bj z(tk)/Bj x(tk)で定義される媒介変数θを用いてsj z(tk) ＝ sinθ、sj x(tk) ＝ cosθと記述される。これを書き直せばsj z(tk) ＝ sin(arctan(Bj z(tk)/ Bj x(tk)))、sj x(tk) ＝ cos(arctan(Bj z(tk)/ Bj x(tk)))である。
　式(9)から明らかなようにBj x(tk)の変数はtkのみでありτとγは定数である。従って、sj z(tk) ＝ sin(arctan(Bj z(tk)/ Bj x(tk)))及びsj x(tk) ＝ cos(arctan(Bj z(tk)/ Bj x(tk)))はBj Z(tk)とtkを変数とする関数としてsj z(tk) ＝ f1(Bj z(tk) , tk)及びsj x(tk) ＝ f2(Bj z(tk) , tk)のような一般化した表現もできる。
　スピンを２次元ベクトルとして記述しているのでsj z(tk)とsj x(tk)の２成分が登場しているが、Bj z(tk)を式(10)に基づき決定するならばsj x(tk)は必要ない。これは、[－1, 1]を値域とするsj z(tk)のみでスピン状態を記述できることに対応している。最終的な解sj zdはsj zd ＝ －1 or 1になる必要があり、sj z(τ) > 0ならばsj zd ＝ 1、sj z(τ) ＜ 0ならばsj zd ＝ －1とする。
　図２に、上記のアルゴリズムをフローチャートにまとめたものを示す。ここでtm ＝τである。図２のフローチャートの各ステップ101~109は、時間t ＝ 0からt ＝ τに到る図1の手順100の、ある時刻での処理に対応している。すなわち、フローチャートのステップ102、104、106がそれぞれ、t ＝ t1, tk＋1, tmにおける上記の式(9)及び(10)に対応している。最終的な解はステップ108において、sj z ＜ 0ならばsj zd ＝ －1、sj z > 0ならばsj zd ＝ 1とすることにより定める（109）。
　ここまでは課題が式(3)で表現された場合に如何に解かれるかを示した。次に具体的課題が如何に局所場gjと変数間相互作用Jij (i, j ＝ 1, 2, …, N)を含む式(3)で表現されるかに関して具体例を挙げて説明する。例えば、具体的課題として、電力供給管理の問題を考える。この場合、局所場は気温のような自然現象の量、あるいは電力使用量とする。すなわち、局所場gj (j ＝ 1 － 10)により各地区の気温を、局所場gj (j ＝ 11 － 20)により各地区にある公共施設（図書館、映画館、スーパーマーケット、等）の電力使用量を、局所場gj (j ＝ 21 － 100)により各民家の電力使用量を表すとする。
　σ^j z (j ＝ 11 － 100)はどこに電力を配分するかを表す変数とする。但し、j ＝ 1 － 10は気温を表す添え字なので、σ^j z (j ＝ 1 － 10)は電力配分を表すのではなく気温を公共施設や民家の活動に影響させるための変数と考える。気温は自然現象で決まるもので人工的要因の影響をほとんど受けないはずなのでσ^j z (j ＝ 1 － 10)が他の変数に影響されないように局所場gj (j ＝ 1 － 10)は大きい値に設定する。
　気温と公共施設及び民家の活動との相関強度は、変数間相互作用Jijを通して表現する。また、気温と電力使用の相関は近年提唱されている電力シェアの概念にも影響される。例えば冷房が必要な時間帯に各家庭ごとにエアコンを使うのではなく公共施設に赴いて各家庭の電力を減らそうとする動きである。その動きは公共施設を表す添え字i ＝ 11 － 20と民家を表す添え字 j ＝ 21 － 100に対して変数間相互作用Jijがゼロでない値を取ることで表現される。但し、この概念に基づく相互作用は気温と民家の活動に関する直接的相関に比べれば小さいものなのでここでの変数間相互作用Jijの値は比較的小さい。また各家庭は独立に生活しているわけではなくお互いに影響を及ぼしあっているはずなので変数間相互作用Jij (i, j ＝ 21 － 100)も有限になる。以上のような考察を通して変数間相互作用Jijを具体的に設定し、式(3)の基底状態探索を通して、最適電力供給配分（σ^j zの固有値 ＝ ＋1あるいは－1）を求める。
　尚、各項目に対するσ^j zを１変数で表現できない場合は複数のσ^j zを使用してよく、それに従い局所場gj及び変数間相互作用Jijも項目ごとに複数使用することになる。σ^j zは電力配分を表す変数であるが、人の動きや公共施設の開館状況と相関している。そのため得られた解によっては「ある公共施設を閉じるべし」との解釈にもなり得る。
　以上が具体的課題を式(3)で表現する簡単な一例である。本実施例が適用できる具体的課題は以上で例示した電力供給管理の問題に限定されず、旅行経路最適化、渋滞回避の為の自動車誘導、回路設計、商品供給管理、スケジューリング、金融資産選択などの多くの課題解決に適用可能である。
　実施例１では量子力学的な式を元に期待値を取ることにより古典量に移行し、図１及び図２を用いて、古典量によるアルゴリズムを説明した。実施例１は本発明のアルゴリズムを説明するのが主な狙いであったために、量子力学的な効果を含めずに説明した。しかし、量子力学的効果を付加すれば正解率向上や演算スピードの向上を期待できる。そこで実施例２では、アルゴリズムそのものは古典的であるが、性能向上のために量子力学に基づく補正パラメタを加える方法を説明する。
　量子力学の特徴としては線形重ね合わせ状態と量子縺れ（非局所相関）がある。例えば、|0>と|1>の２状態を取りうる量子ビットを考える。線形重ね合わせ状態とは|Ψ> ＝ α|0> ＋ β|1>のように２つの状態の和になっているものである。この線形重ね合わせ状態の性質は実施例１でスピンをベクトル的に扱うことによりすでに取り込まれている。即ち、sj z(tk) ＝ 1ならば状態は|0>で、sj z(tk) ＝ －1ならば状態は|1>である。|0>及び|1>はスピンの量子化軸をz軸に選んだ場合の状態に対応し、x軸に向いたsj x(t0) ＝ 1の場合は|Ψ(t0)> ＝ (|0> ＋ |1>)/√2と表現される。またsj x(t) ＝ －1ならば状態は|Ψ(t0)> ＝ (|0> － |1>)/√2である。x軸を考慮したことは線形重ね合わせを考慮したことを意味する。
　本実施例ではもうひとつの量子力学的効果である量子縺れについて述べる。例として２量子ビット系の状態が|Ψ> ＝ α|00> ＋ β|11>と記述できる場合を考える。規格化条件により|α|2 ＋ |β|2 ＝ 1である。|00>及び|11>の１番目の変数が第１量子ビット、２番目の変数が第２量子ビットを表す。パウリのスピン行列の性質としてσ^j z|0> ＝ |0>、σ^j z|1> ＝ －|1>なのでσ^1 z|Ψ> ＝ α|00> － β|11>であり、<Ψ|σ^1 z|Ψ> ＝ |α|2 － |β|2となる。また、σ^1 x|0> ＝ |1>、σ^1 x|1> ＝ |0>なのでσ^1 x|Ψ> ＝ α|10> ＋ β|01>であり<Ψ|σ^1 x|Ψ> ＝ 0となる。また、σ^1 y|0> ＝ i|1>、σ^1 y|1> ＝ －i|0>なのでσ^1 y|Ψ> ＝ iα|10> － iβ|01>であり<Ψ|σ^1 y|Ψ> ＝ 0である。よって<σ^1 x(τ)>2 ＋ <σ^1 y(τ)>2 ＋ <σ^1 z(τ)>2 ＝ (|α|2 － |β|2)2である。極端な例として、量子縺れが最大となるα ＝ βの場合は<σ^1 x(τ)>2 ＋ <σ^1 y(τ)>2 ＋ <σ^1 z(τ)>2 ＝ 0となって第１スピンベクトルの大きさが0になってしまう。このようなことは量子縺れがない場合には起こらない。例えば１スピン系を考え、状態を|Ψ> ＝ α|0> ＋ β|1>とすれば<Ψ|σ^1 z|Ψ> ＝ |α|2 － |β|2、<Ψ|σ^1 x|Ψ> ＝ α*β＋αβ*、<Ψ|σ^1 y|Ψ> ＝ iαβ* － iα*βであり、<σ^1 x(τ)>2 ＋ <σ^1 y(τ)>2 ＋ <σ^1 z(τ)>2 ＝ (|α|2 ＋ |β|2)2 ＝ 1となって確かに大きさが１に保存する。
　以上、一例ではあるが量子縺れがある場合にはスピンベクトルの大きさが１に保存しないことがわかった。古典的にはスピンベクトルの大きさは一定の１のはずだが、量子縺れがあればスピンベクトルは大きさが１にならない。実施例１ではスピンベクトルの大きさが１になることを前提にtanθ ＝ <Bj z(t)>/<Bj x(t)>で定義されるθを媒介変数としてsj z(tk) ＝ sinθ、sj x(tk) ＝ cosθとした。しかし、この方法ではこの系が本来持っている量子縺れの性質を反映したものにはなっていない。そこで量子縺れの性質を反映させることを考える。
　上述のようにスピンベクトルは１に保存しない。そこでスピンの大きさを表す補正パラメタrs (0 ＜＝ rs ＜＝1)を定義する（「＜＝」は「以上」の意）。またスピンベクトルが１に保存しないことに連動して式(8)の比例関係が成り立たなくなる。そこで補正パラメタrBを定義して式(8)を式(11)に変形する。
Figure JPOXMLDOC01-appb-M000011
　実施例１の場合と同様にスピンの向きを表す角度θをtanθ ＝ sj z(tk)/sj x(tk)で定義する。これに式(11)を適用すればtanθ ＝ rB・Bj z(tk)/Bj x(tk)である。スピンの大きさがrsであることを考慮すればsj z(tk) ＝ rs・sinθ、sj x(tk) ＝ rs・cosθとなる。これらの関係式により補正パラメタrs及びrBを通して量子縺れの効果が古典的なアルゴリズムに取り込まれる。θを使わない表記にすればsj z(tk) ＝ rs・sin(arctan(rB・Bj z(tk)/ Bj x(tk)))、sj x(tk) ＝ rs・cos(arctan(rB・Bj z(tk)/ Bj x(tk)))となる。また、rsとrBを関数f1及びf2に組み込めばsj z(tk) ＝ f1(Bj z(tk) , tk)及びsj x(tk) ＝ f2(Bj z(tk) , tk)である。
　この補正パラメタrs及びrBは量子縺れを起源としており、tk, sj z(tk), sj x(tk)に依存して細かく制御されることが好ましい。しかし、量子縺れに関する情報を正確に取得することは原理的に困難であり、何らかの対処法を考える必要がある。実際には問題に応じて半経験的に定めることになるが、大まかな決定法は以下のようになる。rBは符号も変わり得る量で量子縺れを最もよく反映した量である。一方、rsは0 ＜＝ rs ＜＝ 1を満たす補正因子でrBに比べて役割は小さい。従って、全演算時間に亘ってrs ＝~ 1と設定してもよく（「＝~」は「ほぼ等しい」の意）、主としてrBにより量子効果を取り入れる。演算の開始時点では量子縺れは無いのでt ＝ 0においてrB ＝ 1とし、t > 0において徐々にrBを0に近づける。t ＝ τに近づけば多くのスピンはsj z ＝ 1 or －1に収束するが一部のスピンはsj z > 0になるかsj z ＜ 0になるか微妙な振る舞いをする。演算の成否を左右するのはこれら収束性の悪いスピンである。従って、t ＝~ τではこれらのスピンに対して最適になるようにrBを決定する。量子縺れの効果を最大限に取り入れるのでrB ＝~ 0とする。sj z ＝ 1 or －1に収束するスピンは向きが安定しているのでrB ＝~ 0としたことに伴う悪影響はほとんどない。
　以上が時間依存性に関するrBの設定法である。rBに磁場依存性を持たせることも効果的である。Bj z(tk)/Bj x(tk) ＝~ 0の場合は必然的にsj z(tk)/sj x(tk) が不定になる。従って、時刻tの進行に伴いrBがrB ＝~ 1からrB ＝~ 0となる変化をBj z(tk)/Bj x(tk) ＝~ 0の場合には|Bj z(tk)/Bj x(tk)| >> 0の場合に比べて早めることが効果的である。
　また、サイト間の特徴が特別にない場合はrsやrBにサイト依存性を持たせないが、サイト別の特徴が事前にわかっている場合にはそれに対応してrsやrBをサイト依存にすれば解の正解率向上が期待できる。
　図３にrs及びrBを導入した場合のフローチャートを示す。図２のフローチャートとの違いは、ステップ１０３、１０５、１０７がそれぞれ補正パラメタrs及びrBを含むステップ１０３a、１０５a、１０７aに変更された点である。この例では、上述のように、関数fに関して補正パラメタrs及びrBを追加し、tanθ ＝ rB・Bj z(tk)/Bj x(tk)によりθを定義し、sj z(tk) ＝ rs・sinθによってsj z(tk)を定めることとし、従って関数fがf(Bj z(tk) , tk) ＝ rs・sin(arctan(rB・Bj z(tk)/Bj x(tk)))となる。また、補正パラメタrs及びrBは、tkとBj z(tk)に依存して細かく制御されることが望ましい。
　実施例１及び２では課題のハミルトニアンが式(3)で与えられ、各サイトの有効磁場が式(9)及び(10)で与えられるものであった。各サイトには局所場gjがあった。局所場gjがあればsj zはgjで決まる向きが第ゼロ近似の向きになり、Jijで決まる相互作用が加わることによりsj zの向きが補正される。しかし、全サイトがgj ＝ 0となる場合は第ゼロ近似の向きといった概念が存在せず、また一般に縮退数も多くなるためにsj z > 0になるのかsj z ＜ 0になるのか定まらず、演算時間が経過してもsj z ＝~ 0から抜け出せなくなる。また抜け出せたとしても正解でなければスピン間相互作用の結果、スピン間で互いに向きを反転させようとする力が働き、時間のステップごとにスピンの向きが反転する振動現象が生じ、解が収束しない。
　これらの問題を解消するために式(10)に緩和項（ピン止め項）を加えて有効磁場を式(12A)にする。
Figure JPOXMLDOC01-appb-M000012
第３項が緩和項である。sgn(・)は符号関数で、sj z > 0に対してsgn(sj z) ＝ 1、sj z ＝ 0に対してsgn(sj z) ＝ 0、sj z ＜ 0に対してsgn(sj z) ＝ －1である。この緩和項はスピンの向きを元に留めるように働くので上述の振動現象が解消し、解の収束性が向上する。gpinaの値は経験的に決めることになる。緩和項は解の収束性を向上させるための付加的な項であり|Jij|よりも十分に小さい必要がある。一方、小さすぎると十分な働きが期待できない。範囲を規定するとすれば、係数gpinaは|Jij|の平均値の1%から50%の範囲で調整することが適切であると考えられる。ひとつの目安としては、係数gpinaを|Jij|の平均値の1/10程度にすればよい。
　式(12A)の緩和項（第３項）はsj zの符号のみに依存しsj zの大きさには依存しない。一方、第１項はsi zの大きさに依存する。そこで第３項もsj zの大きさに依存させる方法もある。その場合は式(12B)となる。
Figure JPOXMLDOC01-appb-M000013
式(12B)の場合の第３項の大きさはsj zの大きさに依存して変化するので係数gpinbは典型的には|Jij|の平均値程度の大きさであり、範囲で言えば|Jij|の平均値の50%から200%程度になる。
　基底状態が縮退している場合には解のひとつに演算を誘導する必要がある。誘導しなければ解の平均値であるsj z ＝~ 0から抜け出せない。緩和項はこの誘導にも役立つ。解の一つに誘導するためには第ゼロ近似の解を適切に設定する必要がある。そこで基準にするサイト（jサイトとする）を一つ選びsj z(t0) ＝ 1とし、他のサイトの向きをjサイトを基準にJijの符号に基づき初期の段階で定めることにする。このようにすれば適切な第ゼロ近似の設定が為され、その後の局所場応答の演算を通して正解のひとつに収束していく。その際、緩和項は以下のように貢献する。
　t ＝ t0でsj z(t0) ＝ 1, si z(t0) ＝~ 0 (i ≠ j)と設定し、式(12A)あるいは式(12B)と式(11)に基づいて時間発展させる。ここでsi z(t0) ＝~ 0 (i ≠ j)はsi z(t0)が他のサイトにほとんど影響を与えないように例えばsi z(t0) ＝ 1/1000程度の値にすることを意味する。時間ステップt0 → t1により式(12A) あるいは式(12B)と式(11)に基づいてJij ≠ 0となるiサイトに対してはsi z > 0 or si z ＜ 0が得られる。Jij ＝ 0となるiサイトに対してはt0 → t1のステップでsj z(t0) ＝ 1の初期設定は伝播しないが、si z > 0 or si z ＜ 0となったサイトが増えるので次の時間ステップにおいてそのサイトを通してsj z(t0) ＝ 1の初期設定が間接的に伝播し、演算の早い段階でほぼすべてのサイトがsi z > 0 or si z ＜ 0となる。緩和項があるのでjサイトのsj z > 0は演算の初期段階では維持され、jサイト自身にsj z(t0) ＝ 1の履歴が残っている間にほぼすべてのサイトにsj z(t0) ＝ 1の情報が伝播する。これが適切な第ゼロ近似の設定となる。
　図４に緩和項が加わった場合のフローチャートを示す。図４Ａが式(12A)に対応し、図４Ｂが式(12B)に対応する。図３のフローチャートとの違いは、ステップ１０２、１０４、１０６が緩和項を含むステップ１０２b、１０４b、１０６bに変更された点である。
　以上、基底状態が縮退している場合に解のひとつに演算を誘導する方法としてt ＝ t0でsj z(t0) ＝ 1, si z(t0) ＝~ 0 (i ≠ j)と設定することを述べた。この方法では強制的な因子はt ＝ t0におけるsj z(t0) ＝ 1のみである。t > t0では緩和項による履歴のみで強制的因子はない。ひとつの解への誘導性をさらに高めるために全時間で強制的因子を加えるのが有効である。そのためにはひとつのサイトj’に対してだけ付加的に局所場項δgj’を加える方法が考えられる。即ち、式(12A)、(12B)の局所場項gj’にδgj’を加えてgj’→gj’ +δgj’とする。ここで右辺の２つの量は打ち消し合わないようにδgj’の符号はgj’と同符号とする。元々がgj’ = 0ならばδgj’の符号はどちらでも良い。この付加項によりj’サイトのスピンはある向きに強く誘導される。δgj’は|Jij|の平均値の10%から100%程度の範囲で調整することが適切であるが、ひとつの目安としては|Jij|の平均値の50%程度にすればよい。尚、δgj’を通して強制的因子が加わっているので、この場合はsj’ z(t0) ＝ 1の初期設定は不要であり、全サイトに対してsi z(t0) ＝~ 0とすればよい。この場合のフローチャートを図４Ｃ及び図４Ｄに示す。図４Ｃが式(12A)に対応し、図４Ｄが式(12B)に対応する。
　以上説明したように、緩和項を加えたことにより、時間発展におけるスピンの向きに関する振動が抑えられ、解の収束性が向上する。さらに、この緩和項は以下の効果もあった。基底状態の縮退数が多い場合には、系がいずれの基底状態に向かえばよいのか判断できず、解の平均値であるsj z ＝ 0の状態に落ち込み、そこから抜け出せない事態に陥る可能性がある。そこで初期状態においてひとつのサイト（サイトjとする）のスピンだけを明確に定め(sj z ＝ 1 or －1)、その他のスピンをsk z ＝~ 0 (k ≠ j, sk x ＝~ 1)とし、スピン間相互作用を通して時間発展の初期段階でjサイトを基準にしたスピン配列を定める。緩和項があるので演算の初期段階ではjサイトの向きは固定されており、基準が維持される。そのため、縮退解の一つに対応する良好な近似解を演算の初期段階で実現し、そのまま縮退解のひとつに誘導する。このように、緩和項は解の収束性を高めると共に、縮退数が多い場合にも解の一つに収束させる。尚、縮退数が多い場合に解の一つに収束させる方法には２番目の方法もある。式(12A)、(12B)においてひとつのサイトj’に対してのみ付加的に局所場項δgj’を加えてj’サイトのスピンだけをある向きに強く誘導する。これによりj’サイトを基準にした縮退解のひとつに誘導される。
　実施例３では緩和項を導入することにより解の収束性を向上させることを述べた。緩和項が威力を発揮するのは主として全サイトがgj ＝ 0となる場合である。gj ≠ 0となる項があれば緩和項を加えなくても比較的収束性は良い。そこで全サイトがgj ＝ 0であるかどうかに基づき緩和項を加えるかどうかを決定すれば必要な項だけを使った効率的な演算法になる。
　図５に、この場合のアルゴリズムの一例のフローチャートを示す。図５Ａ-図５Ｄはそれぞれ図4Ａ-図４Ｄに対応する。
　上述した式(1)等に見るように、演算時間はτを想定しているが、最終的な解の決定法にはいくつかの方法がある。実施例５として、解の各種決定方法を説明する。
　図６は、最終的な解決定法に関するアルゴリズムの一例をフローチャートとして示す図である。
　第１の方法は図６Ａに示すように、t ＝τ(t ＝ tm)において（115）、sj z > 0ならばsj zd ＝ 1、sj z ＜ 0ならばsj zd ＝ －1とするものである。図２－図５の各実施例のフローチャートはこの場合を記述しており、解決定法のみに焦点を当てたフローチャートが図６Ａである。
　第２の方法は図６Ｂに示すように、sj zの収束性を見るもので、ある時刻tk以降ある十分な時間Δtの間のすべての時刻tqに対してsj zの符号に変化がなければ（121）、その時点のsj zの符号に基づきsj zd ＝ 1 or －1と判定する（122）。
　第３の方法は図６Ｃに示すように、第１の方法と同様にt ＝τ(t ＝ tm)まで演算を継続する（115）が、各段階でのエネルギーを式(3)に基づいて求める。式(3)のσ^j zは固有値が±1 なので演算過程の各段階でsj zの符号に従いσ^j zの固有値が1なのか－1なのかを決める。即ち、sj z(tk) > 0ならばσ^j zの固有値は1 (sj zｄ ＝ 1)、sj z(tk) ＜ 0ならばσ^j zの固有値は－1 (sj z ＝ －1)である。σ^j zの固有値を利用して計算するのはエネルギーに関してであり、演算過程はsj z(tk) (－1 ＜＝ sj z(tk) ＜＝ 1)を利用する。t ＝ τ (t ＝ tm)に到達した時点で各時刻tkでのエネルギーを比較し、最低エネルギーを得た時刻tk’でのsj z(tk’)の符号に基づき最終的な解を決める。
　すなわち、時刻tk各々おいてsj z(tk) ＜ 0ならばsj zd(tk) ＝ －1、sj z(tk) > 0ならばsj zd(tk) ＝ 1として（119）、Hp(tk) ＝ － Σi>jJijsi zd(tk)sj zd(tk) － Σjgjsj zd(tk)を各時刻tkにおいて計算し（123）、Hp(tk)が最小値となった時刻tk’におけるsj zd(tk’)を最終解とする（124）。
　第４の方法は図６Ｄに示すように、第２の方法と同様にすべてのsj zが収束した段階で演算を打ち切る。但し、最終的な解は打ち切った時点でのsj zから判定するのではなく、第３の方法と同様に各時刻でのエネルギーを算出し、最低エネルギーを与えた時刻でのsj z(tk’)から最終的な解を決める（125）。いずれの方法を取るかは使用者が決める。
　図１に示したように時間軸を離散的にした実施例を述べてきた。連続的変化が理想なので時間間隔は細かい程よいが、あまり細かく取りすぎると演算時間が長くなる。そこで時間間隔を演算の進行に合わせて変化させることを考える。
　演算過程で重要な時刻はsj zの符号が変化する時である。演算開始時と終了時近くはsj zの符号変化が比較的少なく、演算の中間段階でsj zの符号の変化が激しい。そこで第１の方法として、プログラム的に演算開始時には時間間隔を大きめにし、時間の経過と共に時間間隔を小さくし、その後逆に時間間隔を大きくしていく設定法がある。
　第２の方法は、各時刻でスピン反転する可能性を評価しその評価に基づき時間間隔を設定する。例えば以下のようにする。すべてのスピンで|sj z|の大きさが概ね等しいならばスピン反転の起こる可能性は低い。この場合は時間間隔を大きくする。一方、特定のスピンの|sj z|の大きさが他のスピンに比べて小さければスピン反転が起こる確率は高い。この場合は時間間隔を小さくする。時間間隔決定法の具体的一例は以下のようにする。最小時間間隔をδtminとする。時刻tkにおける全サイトのスピンの二乗平均をsave(tk)2、最小スピンの二乗の大きさをsmin(tk)2とする。即ち、save(tk)2 ＝ Σj(sj z(tk))2/N, smin(tk)2 ＝ min sj z(tk)2である。[x]をx以下の最大整数としてΔTk＋1,k ＝ tk＋1 － tk ＝δtmin×max(1, [100×(smin(tk)2/save(tk)2)1/2])とする。この例の場合、時間間隔の最小値がδtmin、最大値が100・δtminとなる。
　いずれの方法を取るかは使用者が決めることになる。
　実施例３ではサイトjを選択してsj z(t0) ＝ 1, si z(t0) ＝~ 0 (i ≠ j)とした。サイトjは任意なので、サイトjの選択を変更して同じ問題を解かせることができる。このようにして同じ問題を繰り返し解いて最適解を選択すれば正解率向上を図れる。
　実施例１－７では演算原理及び演算アルゴリズムに関して説明した。実施例８ではまずこのアルゴリズムをプログラムとして動作させる計算機の構成例を説明する。
　図７に本実施例の計算機構成の一例を示す。図７は通常の計算機と基本的に同じ構成で、記憶部である主記憶装置201から、演算部である演算装置202にデータを転送し、演算後、主記憶装置201にデータを返す。これを繰り返して演算を進める。その際の司令塔が制御部としての制御装置203である。本実施例の構成における演算は、図２－図５のフローに従い、時刻別、サイト別に演算装置202で実行する。
　演算装置202で実行されるプログラムは記憶部である主記憶装置201に記憶する。主記憶装置201で記憶容量が足りない場合は、同じく記憶部である補助記憶装置204を利用する。データやプログラム等の入力には入力装置205を使用し、結果の出力には出力装置206を利用する。入力装置205はキーボードのような手入力装置の他、ネットワーク接続のためのインターフェースも含む。また、このインターフェースは出力装置も兼ねる。図７の構成はプログラムとして実施例１－７で説明したアルゴリズムを適用するが、装置そのものは通常の計算機を使う。
　一方、プログラムの実行だけでなく、装置構成も含めて実施例１－７で説明した演算原理、アルゴリズムを利用する方法もある。
　図８にそれを実現した構成図を示す。演算部として局所場応答演算装置1000を含む。これが図７の構成との比較において異なる点である。局所場応答演算装置1000は上述したアルゴリズムを実施するための専用演算部で、実施例1－７で説明した局所場応答演算のみを行い、一般的な処理は演算装置202で実施する。局所場応答演算装置1000で利用する情報は主記憶装置201から転送する。
　必要となる情報は変数間相互作用Jij及び局所場gjといった課題設定パラメタの他、時刻パラメタtkや、量子縺れに係わる補正パラメタrs及びrB、緩和項の係数gpina、gpinb等である。同期を取る等の処理は図７の構成の計算機と同様で制御装置203の役目である。局所場応答演算装置1000からは最終結果だけでなく途中経過の情報が必要に応じて時刻パラメタtkごとに制御装置203の指示に従い、主記憶装置201に転送される。演算の中間段階における局所場応答演算装置1000からの転送値は、例えば演算の中間段階でのイジングスピン・ハミルトニアンのエネルギー算出やsmin(tk)2, save(tk)2の評価に利用する。
　実施例５で説明した第３及び第４の方法では、図６Ｃ、図６Ｄで示したように、中間段階でのイジングスピン・ハミルトニアンのエネルギー値を利用して最終的な解を判定する。sj zdが＋1なのか－1なのかが決まっていればイジングスピン・ハミルトニアンのエネルギー値の計算は単純なものであり、通常の演算装置202を利用する。局所場応答演算装置1000は実施例１－４で述べた局所場応答演算の専用演算部であり、それ以外の処理は通常の演算装置202を利用する。
　実施例８で説明した局所場応答演算装置1000は様々な方法で実現可能である。本実施例ではまず光の並列性を有効利用する方法を述べる。
　図９に全体像を示す。主記憶装置201から演算で必要となるγ, gj, gpina等の情報を制御部1100に、Jijの情報を可変マスク1120に送る。式(12B)対応の場合のgpinbは1120に送り、Jii = gpinbとする。LEDアレイ1110の出力強度は変数sj zの値を反映するものとする。本実施例では光の強度情報のみを利用し、位相情報は利用しない。そのため、光源1110としてLEDアレイのような非干渉性の光源を利用する。干渉性のあるLDを利用することも可能であるが、その場合はLDアレイの出力光同士の干渉効果が現れないように検出器アレイ1130における測定時間を十分に長く取る。LEDアレイ1110からの出力光は横方向のみを広げ、可変マスク1120においてJijに応じて光量を減衰させ、今度は縦方向のみを収束させて検出器アレイ1130で電気信号に変換する。
　図１０は、図９の光学部の一構成例を示す構成図である。LEDアレイ1110から検出器アレイ1130までの光学系は例えば図１０のようなレンズ系を利用する。
　変数sj zは[－1,1]の値を取るが、光源出力は負の値を取りえない。そこでLEDアレイ1110において２つ一組でsj zを表すものとする。即ち、sj z > 0に対してsj z＋ ＝ sj z, sj z－ ＝ 0とし、sj z ＜ 0に対してsj z＋ ＝ 0, sj z－ ＝ sj z ＝ －|sj z|として、LED出力の一方をsj z＋、他方を|sj z－|として、両出力の差sj z ＝ sj z＋ － |sj z－| ＝ sj z＋ ＋ sj z－でsj zを表すものとする。光源側に連動して検出器アレイ1130も２つ一組とする。これにより、やはり負の値を取りえない可変マスク1120に対応できる。検出器アレイ1130で得ようとしている信号はbj z ≡ ΣiJijsi zである。sj zと同様にJij ＝ Jij ＋ ＋ Jij － ＝ Jij ＋ － |Jij －|とすればbj z ＝ ΣiJijsi z ＝ Σi(Jij ＋ ＋ Jij －)(si z＋ ＋ si z－) ＝ Σi(Jij ＋si z＋ ＋ Jij －si z－) ＋ Σi(Jij ＋si z－ ＋ Jij －si z＋)となる。２つ一組とした検出器アレイ1130のそれぞれはbj z＋ ＝ Σi(Jij ＋sj z＋ ＋ Jij －sj z－)及び|bj z－| ＝ Σi(Jij ＋|sj z－| ＋ |Jij －|sj z＋)を検出することになり検出器の差を取ったbj z ＝ bj z＋ － |bj z－| ＝ bj z＋ ＋ bj z－が信号となる。尚、前述のようにJii = gpinbとすれば式(12B)対応となる。
　bj zが得られればgj, gpinaの項を加えることにより式(10)や式(12A)、(12B)に基づきBj zが得られる。この計算は制御部1100で行う。またBj zからsj zを得る計算も制御部1100で行い、sj zの値をLEDアレイ1110に送る。以上で時刻tkからtk＋1への１ステップが終了する。尚、制御部1100は同じ処理を繰り返すものであり、そのための専用回路とする。また、各時刻のsj zは主記憶装置201に転送し解析に利用する。
　図９の実施例ではLEDアレイ1110から検出器アレイ1130までの光学系が自由空間であった。この部分を導波路を使った光学系でも実現可能である。
　図１１にその場合を示す。LEDアレイ（LDアレイ）1110から出力光は分波器1115により分割され、Jijに基づき透過率が設定された可変減衰器1120（図９における可変マスク）を透過後、合波器1125により集光されて検出器アレイ1130で受光される。図９における空間光学系を導波路光学系に変更しただけであり、動作原理は同じである。従って、LEDアレイ（LDアレイ）1110は２つ一組でsj zを表し、検出器アレイ1130も２つ一組で動作する。
　実施例９及び実施例１０では光源を２つ一組にしてsj zを表した。偏光を利用してsj z＋とsj z－を表すことにすればsj zに対して光源をひとつにできる。
　図１２にその場合を示す。光源1111内のLDはそれぞれがsj zを表すためのものである。偏波変調器1112により２つの偏波への分配強度が調整される。sj z＋ ＝ Asin(π/4＋θ), |sj z－| ＝ Acos(π/4＋θ)とすればsj z ＝ sj z＋ ＋ sj z－ ＝ sj z＋ － |sj z－| ＝ √2Asinθでありθ ＜＜ 1とすればsj z ＝~ √2Aθである。θの範囲を－θmax ＜＝ θ ＜＝ θmaxとすればA ＝1/(√2θmax)と選ぶことにより－1 ＜＝ sj z ＜＝ 1が満たされる。偏波変調器1112はsj z＋ ＝ Asin(π/4＋θ), |sj z－| ＝ Acos(π/4＋θ)におけるθを変調する。偏波変調された光は偏光分離器1113により偏波分離され、分波器1115に導かれる。分波器1115から検出器1130までは図１１と同様な動作である。検出器1130は図１１と同様に２つ一組の動作なので差動アンプ1131で差信号を取ってbj z ＝ ΣiJijsi zとする。その後制御部1100に信号を送り実施例９の場合と同様に信号処理がなされ、１ステップ進んだ信号が偏波変調器1112に送られ次のステップに進む。
　本実施例ではθ ＜＜ 1とすることによりsj z ＝ √2Asinθをsj z ＝~ √2Aθとした。ここでθ ＝ arcsinφとすればθ ＜＜ 1の条件を課すことなくsj z ＝ √2Aφとなる。即ち、偏波変調器1112への入力信号の調整により線形性を維持することが可能である。また、実施例２においてt ＝ 0からt ＝ τへの時間変化に伴いrBがrB ＝ 1からrB ＝~ 0に移行することを述べた。この場合にはθ ＝ arcsinφの関数形をさらに変形することになる。
　局所場応答演算装置1000は、実施例９～１１のような光を利用した方法だけでなく電気回路により実現することも可能である。
　図１３にその場合の構成例を示す。各スピンsi zの情報はバッファアレイ1210に一時待機する。メモリ1221にはJijの情報が保持される。式(12B)対応の場合のgpinbは1221で保持し、Jii = gpinbとする。演算器1220ではバッファアレイ1210に待機しているsi zの情報とメモリ1221のJijの情報に基づきbj z ＝ ΣiJijsi zが順次計算され、バッファアレイ1230に転送、保存される。その後制御部1200に信号が転送され、実施例９での制御部1100と同様な信号処理が為されて１ステップ進んだsi zの情報がバッファアレイ1210に送られる。これが１ステップ分の処理であり、この処理がt ＝ t0からt ＝ τまで繰り返される。
　si zは連続量なのでバッファアレイ1210と1230の各si zに対するセルは多ビットで構成し、擬似的な連続量とする。
　本発明における温度の影響は以下のように見積もられる。ビット操作はLED (LD)アレイ1110や偏波変調器1112、バッファアレイ1210及び1230で行われる。ビット反転に必要な電圧は１Vのオーダーである。eを素電荷、kBをボルツマン定数とすれば換算温度はT ＝ eV/kBによりT ＝~ 1.2×104 Kである。この値は室温の300 Kに比べて十分に大きく、実施例９～１２のような構成では温度の影響が無視でき、室温で動作させることができる。
　本発明は上記した実施形態に限定されるものではなく、様々な変形例が含まれる。例えば、ある実施例の構成の一部を他の実施例の構成に置き換えることが可能であり、また、ある実施例の構成に他の実施例の構成を加えることが可能である。また、各実施例の構成の一部について、他の実施例の構成の追加・削除・置換をすることが可能である。
　例えば全数探索を必要とするような課題を処理するための、計算機分野に利用可能である。
100　手順
201　主記憶装置
202　演算装置
203　制御装置
204　補助記憶装置
205　入力装置
206　出力装置
1000　局所場応答演算装置
1100　制御部
1110　LED (LD)アレイ
1111　LD
1112　偏波変調器
1113　偏光分離器
1115　分波器
1120　可変マスク（可変減衰器）
1125　合波器
1130　検出器アレイ
1131　差動増幅器
1200　制御部
1210　バッファアレイ
1220　演算部
1221　メモリ
1230　バッファアレイ

Claims (15)

  1. 　演算部、記憶部、制御部を具備し、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算機であって、
    　N個の変数sj z(j ＝ 1, 2, …, N)が-1 ≦ sj z ≦ 1の値域を取り、局所場gjと変数間相互作用Jij (i, j ＝ 1, 2, …, N)によって課題の設定がなされ、
    　前記演算部では、時刻をm分割して離散的にt ＝ t0 (t0 ＝ 0)からtm (tm ＝ τ)まで演算するものとし、
    　各時刻tkにおける前記変数sj z(tk)を求めるに当たり、前時刻tk-1の変数si z(tk-1) (i ＝ 1, 2, .., N)の値と緩和項の係数gpinaあるいはgpinbを用いてBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ sgn(sj z(tk-1))・gpina}・tk/τあるいはBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ gpinb・sj z(tk-1)}・tk/τを求め、 前記変数sj z(tk)の値域が-1 ≦ sj z(tk) ≦ 1になるように関数fを定めてsj z(tk) ＝ f(Bj z(tk) , tk)とし、
    　時刻ステップをt ＝ t0からt ＝ tmに進めるにつれて前記変数sj zを-1あるいは1に近づけ、
    　最終的にsj z ＜ 0ならばsj zd ＝ -1、sj z ＞ 0ならばsj zd ＝ 1として解を定める、
    ことを特徴とする計算機。
  2. 　前記関数fに関して、
    　ある定数γを用いてBj x(tk) ＝ γ(1 - tk/τ)とし、tanθ ＝ Bj z(tk)/Bj x(tk)によりθを定義し、前記sj z(tk)をsj z(tk) ＝ sinθによって定めることとし、従って前記関数fがf(Bj z(tk) , tk) ＝ sin(arctan(Bj z(tk)/Bj x(tk)))となることを特徴とする請求項１記載の計算機。
  3. 　前記関数fに関して補正パラメタrs及びrBを追加し、
    　tanθ ＝ rB・Bj z(tk)/Bj x(tk)によりθを定義し、sj z(tk) ＝ rs・sinθによって前記sj z(tk)を定めることとし、従って前記関数fがf(Bj z(tk) , tk) ＝ rs・sin(arctan(rB・Bj z(tk)/Bj x(tk)))となることを特徴とする請求項１記載の計算機。
  4. 　前記補正パラメタrs及びrBを、前記tkと前記Bj z(tk)に依存させることを特徴とする請求項３記載の計算機。
  5. 　前記時刻tk各々においてsj z(tk) ＜ 0ならばsj zd(tk) ＝ -1、sj z(tk) ＞ 0ならばsj zd(tk) ＝ 1としてHp(tk) ＝ - Σi＞jJijsi zd(tk)sj zd(tk) - Σjgjsj zd(tk)を各時刻tkにおいて計算し、Hp(tk)が最小値となった時刻tk’におけるsj zd(tk’)を最終解とする、
    ことを特徴とする請求項１記載の計算機。
  6. 　前記係数gpinaは|Jij|の平均値の1%から50%の値であることを特徴とする請求項１記載の計算機。
  7. 　前記課題設定の局所場gj (j ＝ 1, 2, …, N)に関して、
    　gj ≠ 0の成分があれば前記記載のBj z(tk)をBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj}・tk/τにより定め、N個すべてでgj ＝ 0ならば前記記載のBj z(tk)をBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ sgn(sj z(tk-1))・gpina}・tk/τあるいはBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ gpinb・sj z(tk-1)}・tk/τにより定めることを特徴とする請求項１記載の計算機。
  8. 　演算部で実行される演算プログラムであって、
    　N個の変数sj z(j ＝ 1, 2, …, N)が-1 ≦ sj z ≦ 1の値域を取り、局所場gjと変数間相互作用Jij (i, j ＝ 1, 2, …, N)によって課題の設定がなされ、
    　時刻をm分割して離散的にt ＝ t0 (t0 ＝ 0)からtm (tm ＝ τ)まで演算するものとし、
    　各時刻tkにおける変数sj z(tk)を求めるに当たり、前時刻tk-1の変数si z(tk-1) (i ＝ 1, 2, .., N)の値と緩和項の係数gpinaあるいはgpinbを用いてBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ sgn(sj z(tk-1))・gpina}・tk/τあるいはBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ gpinb・sj z(tk-1)}・tk/τを求め、 sj z(tk)の値域が-1 ≦ sj z(tk) ≦ 1になるように関数fを定めてsj z(tk) ＝ f(Bj z(tk) , tk)とし、
    　時刻ステップをt ＝ t0からt ＝ tmに進めるにつれてsj zを-1あるいは1に近づけ、
    　最終的にsj z ＜ 0ならばsj zd ＝ -1、sj z ＞ 0ならばsj zd ＝ 1として解を定める、
    ことを特徴とする演算プログラム。
  9. 　前記関数fに関して、
    　ある定数γを用いてBj x(tk) ＝ γ(1 - tk/τ)とし、tanθ ＝ Bj z(tk)/Bj x(tk)によりθを定義し、前記sj z(tk)をsj z(tk) ＝ sinθによって定めることとし、従って前記関数fがf(Bj z(tk) , tk) ＝ sin(arctan(Bj z(tk)/Bj x(tk)))となることを特徴とする請求８記載の演算プログラム。
  10. 　前記関数fに関して補正パラメタrs及びrBを追加し、
    　tanθ ＝ rB・Bj z(tk)/Bj x(tk)によりθを定義し、sj z(tk) ＝ rs・sinθによって前記sj z(tk)を定めることとし、従って前記関数fがf(Bj z(tk) , tk) ＝ rs・sin(arctan(rB・Bj z(tk)/Bj x(tk)))となることを特徴とする請求項８記載の演算プログラム。
  11. 　前記補正パラメタrs及びrBを、前記tkと前記Bj z(tk)に依存させることを特徴とする請求項１０記載の演算プログラム。
  12. 　前記時刻tk各々おいてsj z(tk) ＜ 0ならばsj zd(tk) ＝ -1、sj z(tk) ＞ 0ならばsj zd(tk) ＝ 1としてHp(tk) ＝ - Σi＞jJijsi zd(tk)sj zd(tk) - Σjgjsj zd(tk)を各時刻tkにおいて計算し、Hp(tk)が最小値となった時刻tk’におけるsj zd(tk’)を最終解とする、
    　ことを特徴とする請求項８記載の演算プログラム。
  13. 　前記係数gpinaは|Jij|の平均値の1%から50%の値であることを特徴とする請求項８記載の演算プログラム。
  14. 　前記課題設定の局所場gj (i, j ＝ 1, 2, …, N)に関して、
    　gj ≠ 0の成分があれば前記記載のBj z(tk)をBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj}・tk/τにより定め、N個すべてでgj ＝ 0ならば前記記載のBj z(tk)をBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ sgn(sj z(tk-1))・gpina}・tk/τあるいはBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ gpinb・sj z(tk-1)}・tk/τにより定めることを特徴とする請求項８記載の演算プログラム。
  15. 　演算部、記憶部、制御部を具備し、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算機であって、
    　当該計算機は、局所場応答演算装置を備え、
    　前記局所場応答演算装置では、
    　N個の変数sj z(j ＝ 1, 2, …, N)が-1 ≦ sj z ≦ 1の値域を取り、局所場gjと変数間相互作用Jij (i, j ＝ 1, 2, …, N)によって課題の設定がなされ、
    　前記演算部では、時刻をm分割して離散的にt ＝ t0 (t0 ＝ 0)からtm (tm ＝ τ)まで演算するものとし、
    　各時刻tkにおける前記変数sj z(tk)を求めるに当たり、前時刻tk-1の変数si z(tk-1) (i ＝ 1, 2, .., N)の値と緩和項の係数gpinaあるいはgpinbを用いてBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ sgn(sj z(tk-1))・gpina}・tk/τあるいはBj z(tk) ＝ {ΣiJijsi z(tk-1) ＋ gj ＋ gpinb・sj z(tk-1)}・tk/τを求め、 sj z(tk)の値域が-1 ≦ sj z(tk) ≦ 1になるように関数fを定めてsj z(tk) ＝ f(Bj z(tk) , tk)とし、
    　時刻ステップをt ＝ t0からt ＝ tmに進めるにつれてsj zを-1あるいは1に近づけ、
    　最終的にsj z ＜ 0ならばsj zd ＝ -1、sj z ＞ 0ならばsj zd ＝ 1として解を定める、
    　処理を行い、
    　当該局所場応答演算装置における処理は、光学系で並列的に処理される複数の光信号の変調を利用して行われるか、もしくは、バッファアレイに蓄積されたデータを並列的に処理することにより行われる、
    ことを特徴とする計算機。
PCT/JP2015/059765 2015-03-27 2015-03-27 計算機、及び演算プログラム WO2016157333A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
US15/561,566 US10318607B2 (en) 2015-03-27 2015-03-27 Computer and computing program
PCT/JP2015/059765 WO2016157333A1 (ja) 2015-03-27 2015-03-27 計算機、及び演算プログラム
JP2017508862A JP6291129B2 (ja) 2015-03-27 2015-03-27 計算機、及び演算プログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2015/059765 WO2016157333A1 (ja) 2015-03-27 2015-03-27 計算機、及び演算プログラム

Publications (1)

Publication Number Publication Date
WO2016157333A1 true WO2016157333A1 (ja) 2016-10-06

Family

ID=57004892

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2015/059765 WO2016157333A1 (ja) 2015-03-27 2015-03-27 計算機、及び演算プログラム

Country Status (3)

Country Link
US (1) US10318607B2 (ja)
JP (1) JP6291129B2 (ja)
WO (1) WO2016157333A1 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2019079225A (ja) * 2017-10-24 2019-05-23 株式会社日立製作所 計算機及び計算方法
WO2020090674A1 (ja) 2018-10-31 2020-05-07 株式会社日立製作所 情報提供装置および情報提供方法
WO2020235300A1 (ja) 2019-05-20 2020-11-26 株式会社日立製作所 ポートフォリオ作成支援装置およびポートフォリオ作成支援方法
WO2021171746A1 (ja) 2020-02-28 2021-09-02 株式会社日立製作所 ポートフォリオ作成支援装置、ポートフォリオ作成支援方法、及びポートフォリオ作成支援システム
US11868967B2 (en) 2019-04-05 2024-01-09 Hitachi, Ltd. Schedule creation assisting device and schedule creation assisting method

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10332023B2 (en) 2017-09-22 2019-06-25 International Business Machines Corporation Hardware-efficient variational quantum eigenvalue solver for quantum computing machines
US10452990B2 (en) * 2017-11-28 2019-10-22 International Business Machines Corporation Cost function deformation in quantum approximate optimization

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009524857A (ja) * 2006-01-27 2009-07-02 ディー−ウェイブ システムズ，インコーポレイテッド 断熱量子計算の方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6445246B2 (ja) * 2014-03-27 2018-12-26 株式会社日立製作所 情報処理装置及び情報処理方法
WO2017037331A1 (en) * 2015-09-04 2017-03-09 Nokia Technologies Oy Method and apparatus for optimization

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009524857A (ja) * 2006-01-27 2009-07-02 ディー−ウェイブ システムズ，インコーポレイテッド 断熱量子計算の方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2019079225A (ja) * 2017-10-24 2019-05-23 株式会社日立製作所 計算機及び計算方法
WO2020090674A1 (ja) 2018-10-31 2020-05-07 株式会社日立製作所 情報提供装置および情報提供方法
US11868967B2 (en) 2019-04-05 2024-01-09 Hitachi, Ltd. Schedule creation assisting device and schedule creation assisting method
WO2020235300A1 (ja) 2019-05-20 2020-11-26 株式会社日立製作所 ポートフォリオ作成支援装置およびポートフォリオ作成支援方法
WO2021171746A1 (ja) 2020-02-28 2021-09-02 株式会社日立製作所 ポートフォリオ作成支援装置、ポートフォリオ作成支援方法、及びポートフォリオ作成支援システム

Also Published As

Publication number Publication date
US20180089144A1 (en) 2018-03-29
JPWO2016157333A1 (ja) 2017-07-13
US10318607B2 (en) 2019-06-11
JP6291129B2 (ja) 2018-03-14

Similar Documents

Publication Publication Date Title
JP6291129B2 (ja) 計算機、及び演算プログラム
JP6483256B2 (ja) 計算機
Wiersema et al. Exploring entanglement and optimization within the hamiltonian variational ansatz
US20200242501A1 (en) Embedding electronic structure in controllable quantum systems
US10789329B2 (en) Systems and methods for removing unwanted interactions in quantum devices
US10846366B1 (en) Selecting parameters for a quantum approximate optimization algorithm (QAOA)
JP7149504B2 (ja) 量子コンピューティング・マシンのための高ハードウェア効率変分量子固有値ソルバを実現するためのシステム、方法、量子コンピューティング・デバイスおよびコンピュータ・プログラム
Salehi et al. Unconstrained binary models of the travelling salesman problem variants for quantum optimization
US20200342345A1 (en) Classification using quantum neural networks
EP3991104A1 (en) Method of computing a solution to a computational problem using a quantum system and apparatus for computing solutions to computational problems
CN116249992A (zh) 囚禁离子量子计算机的优化电路编译器
Fang et al. Minimizing minor embedding energy: an application in quantum annealing
Wang et al. Quantum computing in a statistical context
Lu et al. Optimization of GPU parallel scheme for simulating ultrafast magnetization dynamics model
Faílde et al. Using differential evolution to avoid local minima in variational quantum algorithms
JP6167188B2 (ja) 計算機、及び演算プログラム
CN111527503B (zh) 不均匀量子退火调度
Yu et al. Solving Quantum Many-Particle Models with Graph Attention Network
Mostame et al. Decoherence in a dynamical quantum phase transition
JP6524340B2 (ja) 計算機及び計算方法
Asoudeh et al. Perfect state transfer on spin-1 chains
Robson et al. Advanced quantum Poisson solver in the NISQ era
Wu et al. Finding quantum many-body ground states with artificial neural network
Al‐Rabadi New dimensions in non‐classical neural computing, part II: quantum, nano, and optical
Skotiniotis et al. Quantum resource theory for charge-parity-time inversion

Legal Events

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

Ref document number: 15887486

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2017508862

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 15561566

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 15887486

Country of ref document: EP

Kind code of ref document: A1