JP6524340B2

JP6524340B2 - 計算機及び計算方法

Info

Publication number: JP6524340B2
Application number: JP2018512732A
Authority: JP
Inventors: 辰也戸丸
Original assignee: Hitachi Ltd
Current assignee: Hitachi Ltd
Priority date: 2016-04-22
Filing date: 2016-04-22
Publication date: 2019-06-05
Anticipated expiration: 2036-04-22
Also published as: WO2017183172A1; JPWO2017183172A1

Description

本発明は、全数探索を必要とするような逆問題や組み合わせ最適化問題に対して高速演算を可能にする計算機及び計算機を用いた計算方法に関するものである。

「ビッグデータ」といった言葉に代表されるように現代は巨大なデータ処理が求められている。データ処理法は社会科学的・経済学的な意味も含めて今後発展していくと思われるが、処理そのものは計算機が担う。計算機は初期値を与えてアルゴリズムに基づき計算するものであり順方向の計算は得意であるが、結果から初期値を推定する逆問題や多くの可能性の中から最適解を選ぶ組合せ最適化問題は、最悪の場合に全数探索が必要で一般に苦手である。しかし、様々な処理が必要なビッグデータでは全数探索のような処理も必要不可欠である。

こういった中で近年注目されるようになってきたのが量子アニール、別名断熱量子計算とも呼ばれる手法である（非特許文献１、２）。この方法は、ある物理系の基底状態が解になるように問題を設定し、基底状態を見つけることを通して解を得ようとするものである。問題を設定した物理系のハミルトニアンをH^_pとする。但し、演算開始時のハミルトニアンはH^_pではなく、基底状態に準備しやすい別のハミルトニアンH^₀とする。次に十分に時間を掛けてハミルトニアンをH^₀からH^_pに移行させる。十分に時間を掛ければ系は基底状態に居続け、最終的にハミルトニアンH^_pの基底状態（解状態）を得る。これが量子アニールの原理である。

イジングスピングラスと呼ばれる物理系の基底状態探索法はNP困難と呼ばれる問題にも対応できる（非特許文献３）。また組合せ最適化問題の中で困難度の高い問題はNP困難に属する。さらに計算複雑性理論でPに分類される問題やNPに分類される問題はすべてNP困難問題に帰着できる。よって、イジングスピングラス系で量子アニールを適用すれば組合せ最適化問題をほぼすべて解けることになり、ビッグデータの処理に大きく貢献する。

量子アニールが注目されるもう一つの理由はディコヒーレンスに対して頑強なことである。量子コンピュータでは量子コヒーレンスが計算時間に亘って保たれていなければならなかった。一方、量子アニールでは基底状態が維持されているならば正解が得られる。必ずしも量子コヒーレンスが維持されている必要はない。現状の技術レベルで純粋な量子系を構築することが困難であり、よって量子コヒーレンスを計算時間に亘って維持することが困難であることを考慮すれば量子アニールが注目される理由が理解できる。但し、量子アニールにも欠点がある。量子アニールを実現しうるのは現状では超伝導磁束量子ビット系に限られており（特許文献１、非特許文献４）、極低温冷却装置を必要とするからである。極低温の必要性は実用的なコンピュータ実現のためには課題である。

この課題を解決するために考案された方法が以下で述べる局所場応答法である（特許文献２、非特許文献５）。まず量子アニールを再考する。アニール（焼きなまし）の概念は元々量子・古典に関係なく存在するものであり、量子アニールは量子性を使って古典アニールの性能を向上させようとしたものである。量子アニールにおいて量子コヒーレンスが必ずしも計算時間に亘って維持される必要がなく、基底状態が維持されれば良かったのはそのためである。アニールの概念が量子・古典に関係なく成り立つことを利用すれば量子アニールとは異なる方法論も有り得る。その観点で発明されたのが前述の局所場応答法である。この方法は量子アニールと同様に、演算器としてのスピン系に時刻t = t₀で横磁場を印加し、磁場を徐々に縮小して時刻t ＝ τで解を得るものである。この方法では、演算器そのものは古典的であり、磁場に対するスピンの応答関数に量子力学的情報が付加される。この方法は古典的マシンで動作させるので室温動作が可能であり、極低温が必要な量子アニールの課題を解決するが、量子効果を反映した応答関数を事前に決める必要がある。特許文献２や非特許文献５では経験的、あるいは類似の問題を解いた結果から量子効果を平均的に含む応答関数を決定した。平均的ながら量子効果を含んでいるので純粋に古典的な方法よりも性能が向上した。しかし、性能をさらに向上させるためには個々の問題及び状態に依存した応答関数にする必要がある。

特表2009−524857号公報国際公開WO2015/118639号公報

T.Kadowaki and H.Nishimori, "Quantum annealing in the transverse Ising model," Phys. Rev. E 58, 5355 (1998). E. Farhi, et al., "A quantum adiabatic evolution algorithm applied to random instances of an NP-complete problem," Science 292, 472 (2001). F. Barahona, "On the computational complexity of Ising spin glass models," J. Phys. A: Math. Gen. 15, 3241 (1982). A. P.-Ortiz, "Finding low-energy conformations of lattice protein models by quantum annealing," Scientific Reports 2, 571 (2012). T. Tomaru, "Quasi-Adiabatic Quantum Computing Treated with c-Numbers Using the Local-Field Response," J. Phys. Soc. Jpn. 85, 034802 (2016).

以上述べたように、量子アニールは超伝導磁束量子ビットを用いるために極低温冷却装置を必要とする。また、局所場応答法は室温で動作するものの、性能向上をもたらす量子効果が平均的なものであったために性能向上が限定的なものであった。そこで本発明の目的は、全数探索を必要とするような難しい課題に対して十分な性能持った室温動作可能な計算機を提供することにある。

本発明の一側面は、変数としてのスピンを局所的な有効磁場に応答させる局所場応答法において、時間軸を離散的とし、有効磁場に対するスピンの応答関数を連続する二つの時刻に依存させることにより量子力学的時間発展に類似の時間発展をさせることである。より具体的には以下のようになる。

N個のスピン変数s_j ^z (j = 1, 2, …, N)が−1≦s_j ^z≦1の値域を取り、局所項を表す係数g_jと変数間相互作用を表す係数J_kj (k, j ＝ 1, 2, …, N)によって課題の設定を行い、
時刻をm分割して離散的にt ＝ t₀ (t₀ ＝ 0)からt_m (t_m ＝ τ)まで演算するものとし、
各時刻t_i(i = 1, 2, .., N)では有効磁場に関係する変数B_j ^z0(t_i)、B_j ^z(t_i)とスピン変数s_j ^z(t_i)をこの順番で定めるものとし、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、前記B_j ^z0(t_i)はB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k ^z(t_i−1) + g_j)により定め、B_j ^z(t_i)は0≦u≦ 1を満たすパラメタuを用いてB_j ^z(t_i) = (1−u)B_j ^z0(t_i) + uB_j ^z(t_i−1)とし、B_j ^z(t_i)に因子t_i/τを掛けてB_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τとし、
前記s_j ^z(t_i)は関数fを使ってs_j ^z(t_i) = f(B_eff,j ^z(t_i),t_i)により定めるものとし、該関数fはs_j ^z(t_i)の値域が−1≦s_j ^z(t_i)≦1になるように定義されることを特徴とし、
時刻ステップをt = t₀からt = t_mに進めるにつれて前記変数s_j ^zを−1あるいは1に近づけ、最終的にs_j ^z < 0ならばs_j ^zfd = −1、s_j ^z > 0ならばs_j ^zfd = 1として解を定めることを特徴とする。

本発明の他の一側面は、入力装置、出力装置、記憶装置、一般演算装置、および、局所場応答演算装置を備える計算機としての一面である。ここで、N個の変数s_j ^z (j = 1, 2,…, N)が−1≦s_j ^z≦1の値域を取り、局所項を表す係数g_jと変数間相互作用を表す係数J_kj (k, j ＝ 1, 2, …, N)によって課題の設定を行う。局所場応答演算装置は、時刻をm分割して離散的にt = t₀ (t₀ = 0)からt_m (t_m = τ)まで演算するものとし、各時刻t_i(i ＝ 1, 2, .., N)では変数B_j ^z0(t_i)、B_j ^z(t_i)、s_j ^z(t_i)を順番に定めるものとし、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、B_j ^z0(t_i)はB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k ^z(t_i−1) + g_j)により定め、B_j ^z(t_i)は0≦u≦1を満たすパラメタuを用いてB_j ^z(t_i) = (1-u)Bj^z0(t_i) + uB_j ^z(t_i−1)とし、B_j ^z(t_i)に因子t_i/τを掛けてB_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τとし、s_j ^z(t_i)は関数fを使ってs_j ^z(t_i) = f(B_eff,j ^z(t_i),t_i)により定めるものとし、該関数fはs_j ^z(t_i)の値域が−1≦s_j ^z(t_i)≦1になるように定義される。そして、時刻ステップをt = t₀からt = t_mに進めるにつれて変数s_j ^zを−1あるいは1に近づけ、一般演算装置では、時刻t_i各々においてs_j ^z(t_i) < 0ならばs_j ^zd(t_i) = −1、s_j ^z(t_i) > 0ならばs_j ^zd(t_i) = 1、s_j ^z(t_i) = 0ならばs_j ^zd(t_i) = 0としてH_p(t_k) ＝ − Σ_k>jJ_kjs_k ^zd(t_i)s_j ^zd(t_i) − Σ_jg_js_j ^zd(t_i)を各時刻t_iにおいて計算し、H_p(t_i)が最小値となった時刻t_i’におけるs_j ^zfd = s_j ^zd (t_i’)を最終解とする。

本発明の他の一側面は、演算部、記憶部、制御部を具備し、制御部の制御により、記憶部と演算部との間でデータをやり取りしながら演算を行う計算機を用いた計算方法としての一面である。この方法では、N個の変数s_j ^z (j = 1, 2,…, N)が−1≦s_j ^z≦1の値域を取り、局所項を表す係数g_jと変数間相互作用を表す係数J_kj (k, j ＝ 1, 2, …, N)によって課題の設定を行う。演算部では、時刻をm分割して離散的にt = t₀ (t₀ = 0)からt_m (t_m = τ)まで演算するものとし、各時刻t_i(i ＝ 1, 2, .., N)では変数B_j ^z0(t_i)、B_j ^z(t_i)、s_j ^z(t_i)を順番に定めるものとし、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、B_j ^z0(t_i)はB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k ^z(t_i−1) + g_j)により定め、B_j ^z(t_i)は0≦u≦1を満たすパラメタuを用いてB_j ^z(t_i) = (1-u)B_j ^z0(t_i) + uB_j ^z(t_i−1)とし、B_j ^z(t_i)に因子t_i/τを掛けてB_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τとし、s_j ^z(t_i)は関数fを使ってs_j ^z(t_i) = f(B_eff,j ^z(t_i),t_i)により定めるものとし、該関数fはs_j ^z(t_i)の値域が−1≦s_j ^z(t_i)≦1になるように定義される。そして、時刻ステップをt = t₀からt = t_mに進めるにつれて変数s_j ^zを−1あるいは1に近づけ、最終的にs_j ^z < 0ならばs_j ^zfd = −1、s_j ^z > 0ならばs_j ^zfd = 1として解を定める。

本発明の他の一側面は、上記の計算方法を演算部に実行させるために、記憶部に格納されるソフトウェアたる演算プログラムそのもの、あるいは、これを記憶した記憶媒体としての一面である。

本発明によると、高い正確度で難問を解ける実用的な計算機が実現する。

本発明の原理を模式的に示した概念図である。実施例１に係るアルゴリズムの一例をフローチャートで示した流れ図である。応答関数r_bの具体的な値を示したグラフ図である。実施例２に係るアルゴリズムの一例をフローチャートで示した流れ図である。実施例３に係るアルゴリズムの一例をフローチャートで示した流れ図である。実施例４に係るアルゴリズムの一例をフローチャートで示した流れ図である。実施例５に係るアルゴリズムの一例をフローチャートで示した流れ図である。実施例６に係るアルゴリズムの一例をフローチャートで示した流れ図である。実施例７に係る計算機構成の一例を示したブロック図である。実施例７に係る計算機内の局所場応答演算装置部の詳細に関して一例を示したブロック図である。

以下、演算の原理も交えて、図面に従い本発明の各種実施例を説明する。ただし、本発明は以下に示す実施の形態の記載内容に限定して解釈されるものではない。本発明の思想ないし趣旨から逸脱しない範囲で、その具体的構成を変更し得ることは当業者であれば容易に理解される。

以下に説明する発明の構成において、同一部分又は同様な機能を有する部分には同一の符号を異なる図面間で共通して用い、重複する説明は省略することがある。

本明細書等における「第１」、「第２」、「第３」などの表記は、構成要素を識別するために付するものであり、必ずしも、数または順序を限定するものではない。また、構成要素の識別のための番号は文脈毎に用いられ、一つの文脈で用いた番号が、他の文脈で必ずしも同一の構成を示すとは限らない。また、ある番号で識別された構成要素が、他の番号で識別された構成要素の機能を兼ねることを妨げるものではない。

図面等において示す各構成の位置、大きさ、形状、範囲などは、発明の理解を容易にするため、実際の位置、大きさ、形状、範囲などを表していない場合がある。このため、本発明は、必ずしも、図面等に開示された位置、大きさ、形状、範囲などに限定されない。

実施例１では量子力学的な記述から出発し、それを古典的な形式に移行することを通して本実施例の土台となる原理を述べる。

図１に本実施例の原理を模式的に示す。基本的枠組みは特許文献２及び非特許文献５に記載の局所場応答法と同じである。t = 0において横磁場を印加してスピンを一方向に揃える。その後、横磁場をゆっくりと減少させてt = τで問題設定のハミルトニアンにする。スピンは各時刻でそれぞれに掛かる局所的有効磁場に応答して時間発展する。

問題設定のハミルトニアンとt = 0におけるハミルトニアンをそれぞれ

とし、時刻tにおけるハミルトニアンを

とする。τが演算時間である。１スピン系の類推からサイトjのスピンが受ける有効磁場はB^_eff,j = −∂H^/∂σ^_jで与えられる。

本実施例の局所場応答法は期待値を取った<σ^_j>をスピン変数とみなして古典的マシン上で動作させるものである。式(1)、 (2)から明らかなように<σ^_j>及び<B^_eff,j>はx, z成分だけからなる。そこで応答関数r_b(t)を

のようにx, z成分だけで定義し、スピンの向きをこの応答関数に基づき定める。スピン系が古典的ならば各スピンの応答は各サイトの有効磁場だけで求まり、応答関数はr_b(t)=1になる。しかし、量子力学には非局所相関(entanglement)があり一般にr_b(t)≠1である。すでに言及したように式(5)は期待値をとることにより古典的な式に移行しているが、r_b(t)≠1を通して量子効果が取り込まれる。r_b(t)の値は経験的あるいは類似問題の量子力学的な事前計算により求める。経験的あるいは類似問題を元にするためにここでの量子効果は平均的なものになる。尚、局所場応答法は量子効果を取り入れずにr_b(t)=1として動作させても良い。量子効果を含めなくても局所場応答法自体は動作する。

式(5)の４つの変数<σ^_j ^z(t_i)>、<σ^_j ^x(t_i)>、<B^_eff,j ^z(t_i)>、<B^_eff,j ^x(t_i)>は期待値を取っており古典的な量である。そこで量子力学の記法から古典物理の記法に変更する。即ち、<σ^_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)とする。この記法の変更により式(5)は

となる。

局所場応答法の時間発展は離散的に行い、時刻t_iにおけるB_eff,j ^z(t_i)は式(4)に従い時刻t_i−1におけるs_k ^z(t_i−1)から決める。時刻t_iにおけるs_j ^z(t_i)は式(6)に従い時刻t_iにおけるB_eff,j ^z(t_i)を元に決める。この手続きを繰り返す。図２はこれをフローチャートとしてまとめたものである。

図２の101はアルゴリズムの出発点で初期値の設定を表す。102aでは時刻t_i−1におけるs_k ^z(t_i−1)を使ってB_eff,j ^z(t_i)を求める。また時刻に依存して横磁場強度B_eff,j ^x(t_i)を定める。103ではB_eff,j ^z(t_i)/B_eff,j ^x(t_i)と応答関数r_b(t)によりスピンの向きに相当するtanθ= r_b(t)・B_eff,j ^z(t_i)/B_eff,j ^x(t_i)を求め、スピンの大きさを表すパラメタr_s(t)を使ってs_j ^z(t_i) = r_s(t)・sinθを求める。r_s(t)はr_s(t)²= s_j ^x(t_i)²+ s_j ^z(t_i)²で定義される量で、r_b(t)と同様に経験的あるいは事前計算で決めておく。具体例は２つ後のパラグラフで述べる。図２の102aと103が繰り返し計算の１セットである。このセットをt=τ(=t_m)まで繰り返し、t=τでs_j ^z(t_m)>0ならばs_j ^zfd=1、s_j ^z(t_m)<0ならばs_j ^zfd=−1として解s_j ^zfdを得る（処理201及び202）。

手順103はs_j ^z(t_i) = f(B_eff,j ^z(t_i),t_i)のように関数fを使って一般的に書くこともでき、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)はスピンの大きさを表すパラメタなので0≦r_s(t)≦1である。また、−1≦s_j ^z(t_i)≦1である。ここで、r_b(t)=1及びr_s(t)=1ならば純粋に古典的になる。

図３に応答関数r_b(t)の事前計算の一例を示す。これは８ビット系でJ_ij及びg_jを[−5、5]の一様乱数で決めた場合である。100個の問題の結果であり800点からなる（100個の問題×８ビット）。ここで、B_zx≡<B^_eff,j ^z>/<B^_eff,j ^x>、s_zx≡<σ^_j ^z>/<σ^_j ^x>である。点が厳密に量子力学的に求めた値である。量子力学の非局所相関を反映して応答関数が大きくばらつく。丸は横軸を４０分割して平均を取ったものである。平均化されて応答関数は滑らかなB_zx依存性になる。滑らかならば数個のパラメタで記述可能になる。非特許文献５では４つのパラメタを使って滑らかな応答関数を記述する手法を述べており、図３の実線r_b ⁰(t)はその手法により求めたものである。もうひとつのパラメタr_s(t)も同じ４つのパラメタから求まる。

以上図３により応答関数の一例を示し、非特許文献５を引用してr_b ⁰(t)及びr_s(t)の決定法について言及した。しかし、応答関数r_b ⁰(t)やr_s(t)の決定法はそれに制限されるものではなく様々な方法がありえると共に経験的に定めることもできる。以下の実施例ではr_b ⁰(t)を図３の実線に限定せずに平均的な応答関数を表すものとする。

尚、r_b(t)と同様にr_s(t)もr_s(t)=1として動作させることも可能である。r_b(t)の値域が−∞<r_b(t)<∞であるのに対して0≦r_s(t)≦1なので、r_s(t)=1としたことによる最終解に対する影響はr_b(t)=1としたことによる影響に比べて小さい。従って、r_s(t)決定のための事前情報が不十分な場合にr_s(t)=1に設定することは有効な手段である。

図２のアルゴリズムではt=τ(=t_m)まで繰り返し計算をし、s_j ^z(t_m)の符号に基づき最終的な解s_j ^zfdを定めた。しかし、t=τにおいて最適解が得られるとは限らない。t<τにおいて最適解が得られた後、解の正確度が低下することもある。

図４にその対応としてのアルゴリズムを示す。図４の手順300に示すように、各時刻t=t_iでエネルギーを計算し、全時刻を通して最低エネルギーを与えたスピン配列を最終解にすればよい。

手順300内では各時刻t=t_iでs_j ^z(t_i)の符号判定を行いs_j ^z(t_i)>0ならばs_j ^zd(t_i)=1、s_j ^z(t_i)<0ならばs_j ^zd(t_i)=−1とする。s_j ^z(t_i)=0の場合はs_j ^zd(t_i)=0とする。（処理301及び302）。次にs_j ^zd(t_i)の値に対するエネルギー値H_p(t_i)を求める（手順303）。さらにH_p(t_i)をH_p(t_i−1)と比較する。H_p(t_i)<H_p(t_i−1)ならばt≦t_iの範囲でH_p(t_i)が最低エネルギーなのでt=t_iでのs_j ^zd(t_i)をs_j ^zfdとして保存する。この処理をt=t_mでまで繰り返せばs_j ^zfdに最適解が残る。尚、初期値としてH_p(t₀)=0とする。

図３では応答関数を量子力学的に計算した例を示した。量子力学の非局所相関を反映して応答関数は大きくばらついた。局所場応答法の性能を向上させるためには応答関数にこのばらつきを取り入れる必要がある。本実施例ではその方法を述べる。

局所場応答法の演算において重要な時間帯はスピンの符号が変わる周辺である。その近傍ではs_j ^z(t)≒0であり、それに連動してB_eff,j ^z(t)≒0になる。なお、「≒」は略等しいを意味する。図３における応答関数は式(5)を書き直した

である。<B^_eff,j ^z(t)>が分母にあるので<B^_eff,j ^z(t)>≒0ならばr_b(t)は大きく揺らぐ。局所場応答法は応答関数に従い時間発展させるので、r_b(t)が大きく揺らぐならばそれを再現できるかどうかが解の正確度を左右することになる。

スピンの符号が変化する時刻は量子力学的には線形結合状態になっており、スピンが+1の状態と−1の状態が50:50で重ね合わせなっている。この時刻はバンド理論で言えばバンド反発点である。従って、量子力学的に見てもこの時刻は重要でありs_j ^z(t)≒0、B_eff,j ^z(t)≒0近傍の振る舞いを正しくアルゴリズムに取り入れる必要がある。

線形結合状態は量子力学固有の特徴の一つなので古典的マシンにそのまま導入することはできない。そこで以下に述べるように有効磁場を２時刻のスピン値から決めることにより線形結合と類似の振る舞いをさせることにする。

実施例１では時刻t_iにおける有効磁場B_eff,j ^z(t_i)を時刻t_i−1におけるスピンの値s_j ^z(t_i−1)により求めた。即ち、

が有効磁場B_eff,j ^z(t_i)における基本的因子である。本実施例では時刻t_iだけでなくひとつ前の時刻t_i−1における因子も考慮して

に基づいて有効磁場を決定する。ここでuは0≦u≦1で正確度が高くなるように適当に定める。典型的な値はu≒0.1である。横磁場及びそのスケジュールも含めて有効磁場を記述すれば

である。スピンの符号が反転する時刻周辺（量子力学におけるバンド反発点）ではそれに連動して有効磁場の符号も反転する。そのため式(9)に従う新たな有効磁場ではB_j ^z0(t_i−1)とB_j ^z0(t_i)が相殺してB_j ^z(t_i)≒0になる。スピンの値はこの新たな有効磁場と平均化された応答関数r_b ⁰を用いて

により定める。これによりB_j ^z(t_i)≒0ならばs_j ^z(t)≒0となる。

量子力学ではスピン状態が線形結合状態になっている。一方、ここでの取り扱いでは式(9)で定義した有効磁場が２時刻の有効磁場の線形結合になっている。従って、数式上同じことをしている訳ではない。しかし、式(9)の取り扱いにより量子力学的スピンに類似の振る舞いを古典的マシン上で再現できるようになる。これにより局所場応答法の解の正確度が向上する。

有効磁場を２時刻に依存させたこの取り扱いでは、式(11)に示すように応答関数として平均的な応答関数r_b ⁰を利用した。一方、式(7)に相当する応答関数は、時刻t_iだけで決まる有効磁場B_eff,j ^z0(t_i)=(t_i/τ)B_j ^z0(t_i)を用いて

となる。式(11)と式(12)の定義の違いによりこのr_bはばらつくことになる。式(12)に基づく応答関数は図３に例示したばらついた応答関数を定性的に再現する。

図５に以上の原理に基づくフローチャートを示す。図４との違いは手順102aが手順102bに変わることである。この変更は式(8)、(9)、(10)に基づくものであり有効磁場B_eff,j ^z(t_i)が改良され、応答関数r_b(t_i)のばらついた様子が再現される。また手順103にある応答関数を図４では一般的なr_b(t)により記載したが図５では式(11)に従いr_b ⁰(t)とした。

本実施例ではs_j ^z(t)≒0、B_eff、j ^z(t)≒0近傍の現象を再現するために式(9)の取り扱いを導入した。この取り扱いはこの領域以外でも量子力学的な意味を持つ。簡単な例としてJ₁₂ > 0、 g₁ = g₂ = 0の２スピン系を考える。仮にs₁ ^z(t_i) > 0、 s₂ ^z(t_i) < 0であったとする。このとき式(8)により有効磁場を定めれば相手のスピンだけで自身の有効磁場が定まるので両スピンの符号が反転してs₁ ^z (t_i+1) < 0、 s₂ ^z (t_i+1) > 0となる。以降も同様で両スピンは反転を繰り返す。

しかし、実際の量子力学的スピン系では相手のスピン状態だけでなく自身のスピンの状態にも依存してその後のスピンの状態が定まる。式(9)はこの効果を持つ。式(9)の右辺第１項のB_j ^z0(t_i)は式(8)に基づき時刻t_i−1におけるサイトj以外のスピンの値により定まる。式(9)の右辺第２項のB_j ^z(t_i−1)はs_j ^z(t_i−1)/s_j ^x(t_i−1) = r_bB_j ^z(t_i−1)/B_j ^x(t_i−1)の関係を通して時刻t_i−1におけるサイトj（自身のサイト）のスピンの値に依存している。従って、式(9)に基づく有効磁場決定法はt = t_i−1での相互作用相手及び自身の両者に依存した形になり、両者に依存している点で量子力学的である。具体的な効果を上記の２スピン系の例で見れば、逆符号となるB_j ^z0(t_i)とB_j ^z0(t_i−1)が上手く相殺されて不自然な振動が解消される

量子力学的には有効磁場は式(4)に基づき定まる。σ^_k ^zの固有値は±1である。しかし局所場応答法ではスピン変数s_k ^zが期待値<σ^_k ^z>の値を取るように動作させるので|s_k ^z|≦1である。従って、一般にg_jに比べてΣ_k(≠j)J_kjs_k ^zの項を過小評価することになる。実施例３の式(9)における取り扱いでもこの事情に変わりはない。Σ_k(≠j)J_kjs_k ^zの項を過小評価したまま演算させると解の正確度が下がる。そこでs_k ^zの値を参考にg_jの値を規格化することにする。g_jに因子c_i=(Σ_ks_k ^z(t_i−1)²/N)^1/2を掛けてg_j ^norm(t_i) = c_ig_jとすればg_j ^norm(t_i)とΣ_k(≠j)J_kjs_k ^zの項の寄与が概ね同等になり、解の正確度が向上する。但し、離散的に扱う時間軸の分割数をm（従って、t_m=τ）としてc₁=1/mとする。これはc_i=(Σ_ks_k ^z(t_i−1)²/N)^1/2に基づけばs_k ^z(t₀)=0によりc₁=0となってしまうことに対処するためである。

図６に以上の取り扱いを含めたフローチャートを示す。図５との違いは手順102bが手順102cに変わることである。手順102cでは因子c_n=(Σ_ks_k ^z(t_i−1)²/N)^1/2の取り扱いが加わっている。

実際の量子力学的スピン系では常にスピン同士で影響し合っている。即ち、あるサイトjのスピンσ^_j ^zは別のサイトkのスピンσ^_k ^zに影響し、逆にσ^_k ^zがσ^_j ^zに影響する。従って、スピンσ^_j ^zはサイトkのスピンσ^_k ^zを経由して自身に影響する。量子力学においてあるスピンの状態が相互作用相手のスピンの状態だけでなく自身のスピンの状態に依存するのはそのためである。相互作用を通した自身への影響の大きさはΣ_k(≠j)J_kj ²に比例する。ところで応答関数がr_b(t)≠1となったのは相互作用があったことによる。相互作用が無ければr_b(t)=1である。そこで、δr_b(t)≡1−r_b(t)∝Σ_k(≠j)J_kj ²と考えることにする。これにより量子効果が平均的であった応答関数r_b ⁰(t)を個別の問題に依存させることができる。具体的には1−r_b ⁰ _mod(t)=(1−r_b ⁰(t))Σ_k(≠j)J_kj ²/Σ_k(≠j)ave(J_kj ²)で定義される応答関数r_b ⁰ _mod(t)を平均的応答関数r_b ⁰(t)に換えて利用する。ここでave(J_kj ²)はr_b ⁰(t)を決定する際に使用した問題のJ_kj ²を平均したものである。この改良により応答関数が実際の問題に最適化されるようになり解の正確度が向上する。

図７に以上の取り扱いを含めたフローチャートを示す。図６との違いは手順103が手順103cに変わることである。

尚、相互作用を通した自身への影響の大きさがΣ_k(≠j)J_kj ²に比例することは３次の摂動論を用いて理論的に示すことができる。３次の摂動論において有効磁場とスピンの期待値は時刻パラメタη≡t/τを用いて以下になる。

添え字(0−3)は0 − 3次の摂動項を含むことを表す。両式の右辺第１項から第３項までは同じであるが式(14)には第４項が加わる。第４項は第３項でk = jとしたものに−1/2を掛けたもので、この項がサイトiを媒介にして自身のサイトjに影響する項である。この項の符号が負である点が重要で、これにより第１項が部分的にキャンセルされる。第２項と第３項は平均的にはゼロなので平均的振る舞いとしてr_b < 1になる。

式(14)の第４項はΣ_i(≠j)J_ij ²に比例している。即ち、1−r_b(t)∝Σ_i(≠j)J_ij ²である。

実施例５まで、量子力学的考察を通して解の正確度を向上させる方法を述べたが、これらの方法で必ず正解に辿り着くとは限らない。そこで本実施例では新たに補助的な手段を導入する。本実施例の方法では横磁場強度を減少することを通して演算する。本実施例ではこの過程をn_r回繰り返す。

図８にそのためのアルゴリズムの一例を示す。横磁場の減少印加の１ループが手順10である。それぞれのループ内では102cと103cの手順を行い、i→i＋1としてそれを繰り返す。t_i=t_thnになればループを抜け出しスピンの初期値を設定し直して次のループに進む。スピンの初期値は１回目のループ(n=1)ではs_j ^z(t₀)=0とする。２回目以降のループはそれまでに得られたスピンの値を参考にして決める。

例えば２回目のループの初期値は、１回目のループ終了時の最適解s_j ^zfdを使ってs_j ^z(t_i)=s_j,th1 ^z= s_j ^zfd/mとする（手順111）。３回目のループの初期値は例えば、２回目のループ終了時の最適解s_j ^zfdを使うと共にスピンの符号を反転する。即ち、s_j ^z(t_i)=s_j,th2 ^z=−s_j ^zfd/mとする。４回目のループの初期値は例えば、３回目のループ終了時の値を２値化して利用する。即ち、s_j ^z(t_i)=s_j,th3 ^z= s_j ^zfd(t_i−1)/m。５回目のループの初期値は例えば、２回目のループの場合と同様に４回目のループ終了時の最適解s_j ^zfdを使ってs_j ^z(t_i)=s_j,th4 ^z= s_j ^zfd/mとする。このように色々な方法で初期値を変更すれば解の探索空間が広がり正解に辿り着く確率が向上する。５回目のループで演算を終了するならばn_r=5である。

手順102cに示すように有効磁場の強度は

である。横磁場の減少印加をn_r回実施するために１回のループの時間範囲がτではなくt_thn - t_th(n−1)になる。係数a_nはそのために導入した係数で横磁場強度の変化量を調整する。a_n(t_thn - t_th(n−1)) = τと設定すれば100%横磁場を振ることになる。但し、必ずしも100%振る必要はなく、例えばa_n(t_thn - t_th(n−1)) = 0.6τに設定すれば横磁場強度の変化量を60%に抑えることになる。

上述の例では３回目のループでスピン反転を行った。これは局所最適解に落ち込んだスピン状態をスピン配列空間上で最も遠い配置に設定し、大域最適解を探し易くするためである。上述の例では３回目のループの初期値としてスピン反転を行ったが、ループ内の途中でスピン反転することも有効である。この場合、時刻t_i = t_invnを定義しておき、スピン反転のタイミングを制御すればよい。

実施例４においてg_j(t_i)の項とΣ_k(≠j)J_kjs_k ^zの項を同等に寄与させるためにg_j(t_i)をg_j ^norm(t_i) = c_ig_jに置き換えることを述べた。この処理はスピンの初期値がs_j ^z(t₀)=0の場合に最適である。しかし、本実施例における多ループ構成の場合には２回目以降のループで、そこまでの最適解s_j ^zfdを利用して初期値を定めた。この場合、初期値のs_k ^zが最適解候補なのでs_j ^z(t₀)=0の場合に比べてΣ_k(≠j)J_kjs_k ^zに対する過小評価の度合いが縮小する。そのため、２回目以降のループではg_j ^norm(t_i)に因子f(≧1)を掛けてfg_j ^norm(t_i)とし、g_j ^norm(t_i)の寄与を拡大することが有効である。fの値は2から50程度で経験的に定める。

本実施例はアルゴリズムとして示されており、通常のコンピュータ上でソフトウェアとして動作させることも専用のハードウェア上で動作させることもできる。本実施例の特徴は演算が比較的単純であることと並列性が高いことが挙げられる。従って、専用ハードウェアを構築する場合は並列性の高い構成にすることが好ましい。既存のハードウェアを利用するならばGPGPU (General-Purpose computing on Graphics Processing Units)のような並列性の高いハードウェアを利用するのが効果的である。本実施例では本発明を効果的に動作させるための装置の構成例を示す。

図９に本実施例の計算機構成の一例を示す。図９は通常の計算機の構成と類似であるが局所場応答演算装置600を含む。局所場応答演算装置600は実施例１−６で述べた演算を専門に行う部分であり、その他の一般的演算は一般演算装置502で行う。

以上の構成は、単体のコンピュータで構成してもよいし、あるいは、主記憶装置501、一般演算装置502、制御装置503、補助記憶装置504、入力装置505、出力装置506等の任意の部分がネットワークで接続された他のコンピュータで構成してもよい。

一般的な演算は通常の計算機と同様な手順で動作させる。記憶部である主記憶装置501と演算部である一般演算装置502間でデータをやり取りし、その繰り返しで演算を進める。その際の司令塔が制御部としての制御装置503である。一般演算装置502で実行されるプログラムは記憶部である主記憶装置501に記憶させる。主記憶装置501で記憶容量が足りない場合は、同じく記憶部である補助記憶装置504を利用する。データやプログラム等の入力には入力装置505を使用し、結果の出力には出力装置506を利用する。入力装置505はキーボードのような手入力装置の他、ネットワーク接続のためのインターフェースも含む。また、このインターフェースは出力装置も兼ねる。

本実施例の局所場応答演算では、実施例１−６で述べたようにN個のスピン変数s_j ^z(t)とN個の有効磁場変数B_eff,j ^z(t)を交互に繰り返し求める。この繰り返し演算を専門的に実施するのが局所場応答演算装置600である。一方、実施例２で述べた各時刻のエネルギー計算等、繰り返し演算以外の演算は一般演算装置502で行う。

図１０は局所場応答演算装置部の詳細を記したのものである。局所場応答演算装置600は専用記憶装置601、演算部611、612、613、レジスタ621、622、623からなる。演算部611、612、613は互いに独立に動作する並列演算装置であり、それぞれがサイト別に有効磁場B_eff,j ^z(t)とスピン変数s_j ^z(t)を計算する。計算結果のs_j ^z(t)は専用記憶装置601に保存される。演算部611、612、613で必要になる情報J_ij、g_jやs₁ ^z(t)、s₂ ^z(t)、…、s_N ^z(t)はレジスタ621、622、623から読み出す。レジスタ621、622、623はそれぞれ演算部611、612、613に直結した専用記憶部であり、これらレジスタがあるおかげで演算部611、612、613がそれぞれ独立に高速処理できる。図10では各レジスタを便宜上a領域とb領域に分類してある。演算を通して値の変わらないJ_ij、g_jを記憶しているのがa領域で、演算に伴い値が変化するs₁ ^z(t)、s₂ ^z(t)、…、s_N ^z(t)を記憶しているのがb領域である。a領域に保存されるJ_ij、g_jは予め主記憶装置501から専用記憶装置601へ転送し、さらにレジスタ621a、622a、623aに転送しておく。b領域の621b、622b、623bに記憶されるs₁ ^z(t)、s₂ ^z(t)、…、s_N ^z(t)は専用記憶装置601から転送される。このデータ転送はs₁ ^z(t)、s₂ ^z(t)、…、s_N ^z(t)をまとめて送るだけでありランダムアクセスの必要性がないため高速である。以上、レジスタ621、622、623と演算部611、612、613の構成により並列性と高速性が実現される。

実施例2で述べたように、演算では有効磁場B_eff,j ^z(t)とスピン変数s_j ^z(t)を繰り返し求めるだけでなく各時刻でs_j ^z(t)を2値化し、それを元にエネルギーを計算する。この計算は専用記憶装置601から主記憶装置501にデータ転送して一般演算装置502を利用する。即ち、s_j ^z(t)とB_eff,j ^z(t)の繰り返し計算以外は一般演算装置502を利用する。これにより最も計算時間を要する繰り返し計算が効率的になる。

組合せ最適化問題の中で困難度の高い問題はNP困難に属する。またPに分類される問題やNPに分類される問題はすべてNP困難問題に帰着できる。よって、NP困難な組合せ最適化問題を解ければほぼすべての組合せ最適化問題を解けることになる。式(1)の基底状態探索問題はNP困難問題にも対応可能である。本実施例ではその対応の様子を代表的なNP困難問題である最大カット問題を例にして示す。

最大カット問題とはグラフ理論の問題である。グラフ理論ではグラフGを頂点集合Vと辺集合Eから構成しG = (V, E)と書く。辺eは２つの頂点を利用してe = {i, j}と書く。辺eに向きを含めて定義するグラフを有向グラフ、向きの定義を含めないグラフを無向グラフと言う。辺eには重みも定義されておりそれをw_ij、 w_jiと書く。無向グラフならばw_ij = w_jiである。MAX−CUT問題とは、重み付き無向グラフG = (V, E)の頂点を２つのグループに分ける問題において、カットされる辺の重みの総和を最大化する分割法を求める問題である。分割後の２つの無向グラフをG₁ = (V₁, E₁)及びG₂ = (V₂, E₂)とすればMAX−CUT問題は

を最大化する問題である。頂点i∈V₁に対してs_i = 1、頂点j∈V₂に対してs_j = −1とすれば

と書ける。最右辺第１項はグラフGが定まれば定数なのでMAX−CUT問題はΣ_i>js_is_jを最小化する問題となる。イジングスピングラスのハミルトニアンは

で与えられる。よってMAX−CUT問題はJ_ij = −w_ij、 g_j = 0とした式(1)の基底状態探索問題と等価になる。

以上のべた実施例は、は古典的マシン上で動作させるものであり、極低温にする必要が無く、また量子コヒーレンスを考慮する必要もない。その結果、使用可能なリソースが広範囲になり電気回路等も利用できる。さらに、応答関数を連続する二つの時刻に依存させたことにより解の正確度が向上する。これらの性質により高い正確度で難問を解ける実用的な計算機が実現する。

本発明は上記した実施形態に限定されるものではなく、様々な変形例が含まれる。例えば、ある実施例の構成の一部を他の実施例の構成に置き換えることが可能であり、また、ある実施例の構成に他の実施例の構成を加えることが可能である。また、各実施例の構成の一部を、他の実施例の構成に追加・置換することや他の実施例の構成から削除することが可能である。

各種データの解析処理に応用が可能である。

10 〜 300 フローチャートにおける各手順を表す
501 主記憶装置
502 一般演算装置
503 制御装置
504 補助記憶装置
505 入力装置
506 出力装置
600 局所場応答演算装置
601 専用記憶装置
611, 612, 613 演算部
621a, 621b, 622a, 622b, 623a, 623b レジスタ

Claims

演算部、記憶部、制御部を具備し、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算機であって、
N個の変数s_j ^z (j = 1, 2,…, N)が−1≦s_j ^z≦1の値域を取り、局所項を表す係数g_jと変数間相互作用を表す係数J_kj (k, j ＝ 1, 2, …, N)によって課題の設定を行い、
前記演算部では、時刻をm分割して離散的にt = t₀ (t₀ = 0)からt_m (t_m = τ)まで演算するものとし、
各時刻t_i(i ＝ 1, 2, .., N)では変数B_j ^z0(t_i)、B_j ^z(t_i)、s_j ^z(t_i)を順番に定めるものとし、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、前記B_j ^z0(t_i)はB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k ^z(t_i−1) + g_j)により定め、前記B_j ^z(t_i)は0≦u≦1を満たすパラメタuを用いてB_j ^z(t_i) = (1-u)B_j ^z0(t_i) + uB_j ^z(t_i−1)とし、B_j ^z(t_i)に因子t_i/τを掛けてB_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τとし、
前記s_j ^z(t_i)は関数fを使ってs_j ^z(t_i) = f(B_eff,j ^z(t_i),t_i)により定めるものとし、該関数fはs_j ^z(t_i)の値域が−1≦s_j ^z(t_i)≦1になるように定義され、
時刻ステップをt = t₀からt = t_mに進めるにつれて前記変数s_j ^zを−1あるいは1に近づけ、最終的にs_j ^z < 0ならばs_j ^zfd = −1、s_j ^z > 0ならばs_j ^zfd = 1として解を定めることを特徴とする計算機。
前記関数fに関して、
ある定数γを用いてB_eff,j ^x(t_i) = γ(1 − t_i/τ)としてtanθ = B_eff,j ^z(t_i)/B_eff,j ^x(t_i)によりθを定義し、前記s_j ^z(t_i)をs_j ^z(t_i) = sinθによって定めることとし、従って前記関数fがf(B_eff,j ^z(t_i),t_i) = sin{arctan(B_eff,j ^z(t_k)/B_eff,j ^x(t_k))}となることを特徴とする請求項１記載の計算機。
前記関数fに関して補正パラメタr_s及びr_bを追加し、
tanθ = r_b・B_eff,j ^z(t_i)/B_eff,j ^x(t_i)によりθを定義し、s_j ^z(t_i) = r_s・sinθによって前記s_j ^z(t_i)を定めることとし、従って前記関数fがf(B_eff,j ^z(t_i), t_i) ＝ r_s・sin{arctan(r_b・B_eff,j ^z(t_k)/B_eff,j ^x(t_k))}となることを特徴とする請求項２記載の計算機。
前記時刻t_i各々においてs_j ^z(t_i) < 0ならばs_j ^zd(t_i) = −1、s_j ^z(t_i) > 0ならばs_j ^zd(t_i) = 1、s_j ^z(t_i) = 0ならばs_j ^zd(t_i) = 0としてH_p(t_k) ＝ − Σ_k>jJ_kjs_k ^zd(t_i)s_j ^zd(t_i) − Σ_jg_js_j ^zd(t_i)を各時刻t_iにおいて計算し、H_p(t_i)が最小値となった時刻t_i’におけるs_j ^zfd = s_j ^zd (t_i’)を最終解とすることを特徴とする請求項１記載の計算機。
前記g_jを用いてg_j ^norm(t_i) = (Σ_k(s_k ^z(t_i-1))²/N)^1/2・g_jを定義し、前記B_j ^z0(t_i)をB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k(t_i−1) + g_j ^norm(t_i))により定めることを特徴とする請求項１記載の計算機。
前記補正パラメタr_bに対してδr_b≡1-r_bを定義し、δr_b(t)∝Σ_k(≠j)J_kj ²とすることを特徴とする請求項３記載の計算機。
n_r+1個の時刻t_thn (n = 0, 1, 2, …, n_r)とn_r個の係数a_n (n = 1, 2, …, n_r)を定め、t_th0 = 0及びt_thnr =τとし、前記B_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τをt_th(n−1)≦t < t_thnの時間領域においてB_eff,j ^z(t_i) = B_j ^z(t_i)a_n(t_i − t_thn)/τとすることを特徴とする請求項１記載の計算機。
前記時刻t_i各々においてs_j ^z(t_i) < 0ならばs_j ^zd(t_i) = −1、s_j ^z(t_i) > 0ならばs_j ^zd(t_i) = 1、s_j ^z(t_i) = 0ならばs_j ^zd(t_i) = 0としてH_p(t_k) ＝ − Σ_k>jJ_kjs_k ^zd(t_i)s_j ^zd(t_i) − Σ_jg_js_j ^zd(t_i)を各時刻t_iにおいて計算し、H_p(t_i)< H_p(t_i−1)ならばs_j ^zfd = s_j ^zd (t_i)とし、t < t_thnの時間領域で得られたスピンの値s_j ^zfdを利用して前記時刻t = t_thnでのスピンの値s_j ^z(t_thn)を定めることを特徴とする請求項７記載の計算機。
n_q個の時刻t_invn (n = 1, 2, …, n_q)を定め、時刻t_i = t_invnではB_j ^z0(t_i) = (−Σ_k(≠j)J_kjs_k(t_i−1) + g_j)及びB_j ^z(t_i) = (1−u)B_j ^z0(t_i) − uB_j ^z(t_i−1)とすることを特徴とする請求項１記載の計算機。
前記補正パラメタr_s及びr_bを、前記t_iと前記B_j ^z(t_i)に依存させることを特徴とする請求項３記載の計算機。
入力装置、出力装置、記憶装置、一般演算装置、および、局所場応答演算装置を備え、
N個の変数s_j ^z (j = 1, 2,…, N)が−1≦s_j ^z≦1の値域を取り、局所項を表す係数g_jと変数間相互作用を表す係数J_kj (k, j ＝ 1, 2, …, N)によって課題の設定を行い、
前記局所場応答演算装置では、
時刻をm分割して離散的にt = t₀ (t₀ = 0)からt_m (t_m = τ)まで演算するものとし、
各時刻t_i(i ＝ 1, 2, .., N)では変数B_j ^z0(t_i)、B_j ^z(t_i)、s_j ^z(t_i)を順番に定めるものとし、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、前記B_j ^z0(t_i)はB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k ^z(t_i−1) + g_j)により定め、前記B_j ^z(t_i)は0≦u≦1を満たすパラメタuを用いてB_j ^z(t_i) = (1-u)B_j ^z0(t_i) + uB_j ^z(t_i−1)とし、B_j ^z(t_i)に因子t_i/τを掛けてB_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τとし、
前記s_j ^z(t_i)は関数fを使ってs_j ^z(t_i) = f(B_eff,j ^z(t_i),t_i)により定めるものとし、該関数fはs_j ^z(t_i)の値域が−1≦s_j ^z(t_i)≦1になるように定義され、
時刻ステップをt = t₀からt = t_mに進めるにつれて前記変数s_j ^zを−1あるいは1に近づけ、
前記一般演算装置では、
前記時刻t_i各々においてs_j ^z(t_i) < 0ならばs_j ^zd(t_i) = −1、s_j ^z(t_i) > 0ならばs_j ^zd(t_i) = 1、s_j ^z(t_i) = 0ならばs_j ^zd(t_i) = 0としてH_p(t_k) ＝ − Σ_k>jJ_kjs_k ^zd(t_i)s_j ^zd(t_i) − Σ_jg_js_j ^zd(t_i)を各時刻t_iにおいて計算し、
H_p(t_i)が最小値となった時刻t_i’におけるs_j ^zfd = s_j ^zd (t_i’)を最終解とすることを特徴とする計算機。

計算機。
演算部、記憶部、制御部を具備し、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算機を用いた計算方法であって、
N個の変数s_j ^z (j = 1, 2,…, N)が−1≦s_j ^z≦1の値域を取り、局所項を表す係数g_jと変数間相互作用を表す係数J_kj (k, j ＝ 1, 2, …, N)によって課題の設定を行い、
前記演算部では、時刻をm分割して離散的にt = t₀ (t₀ = 0)からt_m (t_m = τ)まで演算するものとし、
各時刻t_i(i ＝ 1, 2, .., N)では変数B_j ^z0(t_i)、B_j ^z(t_i)、s_j ^z(t_i)を順番に定めるものとし、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、前記B_j ^z0(t_i)はB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k ^z(t_i−1) + g_j)により定め、前記B_j ^z(t_i)は0≦u≦1を満たすパラメタuを用いてB_j ^z(t_i) = (1-u)B_j ^z0(t_i) + uB_j ^z(t_i−1)とし、B_j ^z(t_i)に因子t_i/τを掛けてB_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τとし、
前記s_j ^z(t_i)は関数fを使ってs_j ^z(t_i) = f(B_eff,j ^z(t_i),t_i)により定めるものとし、該関数fはs_j ^z(t_i)の値域が−1≦s_j ^z(t_i)≦1になるように定義され、
時刻ステップをt = t₀からt = t_mに進めるにつれて前記変数s_j ^zを−1あるいは1に近づけ、最終的にs_j ^z < 0ならばs_j ^zfd = −1、s_j ^z > 0ならばs_j ^zfd = 1として解を定めることを特徴とする計算方法。
前記g_jを用いてg_j ^norm(t_i) = (Σ_k(s_k ^z(t_i-1))²/N)^1/2・g_jを定義し、前記B_j ^z0(t_i)をB_j ^z0(t_i) = (Σ_k(≠j)J_kjs_k(t_i−1) + g_j ^norm(t_i))により定めることを特徴とする請求項１２記載の計算方法。
n_r+1個の時刻t_thn (n = 0, 1, 2, …, n_r)とn_r個の係数a_n (n = 1, 2, …, n_r)を定め、t_th0 = 0及びt_thnr =τとし、前記B_eff,j ^z(t_i) = B_j ^z(t_i)・t_i/τをt_th(n−1)≦t < t_thnの時間領域においてB_eff,j ^z(t_i) = B_j ^z(t_i)a_n(t_i − t_thn)/τとすることを特徴とする請求項１２記載の計算方法。
前記関数fに関して、
ある定数γを用いてB_eff,j ^x(t_i) = γ(1 − t_i/τ)としてtanθ = B_eff,j ^z(t_i)/B_eff,j ^x(t_i)によりθを定義し、前記s_j ^z(t_i)をs_j ^z(t_i) = sinθによって定めることとし、従って前記関数fがf(B_eff,j ^z(t_i),t_i) = sin{arctan(B_eff,j ^z(t_k)/B_eff,j ^x(t_k))}となり、
前記関数fに関して補正パラメタr_s及びr_bを追加し、
tanθ = r_b・B_eff,j ^z(t_i)/B_eff,j ^x(t_i)によりθを定義し、s_j ^z(t_i) = r_s・sinθによって前記s_j ^z(t_i)を定めることとし、従って前記関数fがf(B_eff,j ^z(t_i), t_i) ＝ r_s・sin{arctan(r_b・B_eff,j ^z(t_k)/B_eff,j ^x(t_k))}となり、
前記補正パラメタr_bに対してδr_b≡1-r_bを定義し、δr_b(t)∝Σ_k(≠j)J_kj ²とすることを特徴とする請求項１２記載の計算方法。