JP6967422B2 - 計算機及び計算方法 - Google Patents

Description

本発明は、全数探索を必要とするような組み合わせ最適化問題等に対して高速演算を可能にする計算技術に関するものである。

IoT (Internet of Things)といった言葉に代表されるように、現代はあらゆるモノがインターネットにつながり，モノから情報が収集され，総合した情報を元にモノが制御される。制御では多くの選択肢の中から最適解が見出されて実行される。極論すれば現代の情報技術は最適解を見出すものとも言える。

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

イジング・スピングラスと呼ばれる物理系を利用した基底状態探索法はNP困難と呼ばれる問題にも対応できる。また組合せ最適化問題の中で困難度の高い問題はNP困難に属する。さらに計算複雑性理論でPに分類される問題やNPに分類される問題はすべてNP困難問題に帰着できる。よって、イジング・スピングラス系で量子アニールを適用すれば組合せ最適化問題をほぼすべて解けることになり、情報技術の最も重要な課題が解決される。

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

この課題を解決するために考案された方法が以下で述べる局所場応答法である（特許文献１〜３、非特許文献１、２）。まず量子アニールを再考する。アニール（焼きなまし）の概念は元々量子・古典に関係なく存在するものであり、量子アニールは量子性を使って古典アニールの性能を向上させようとしたものである。量子アニールにおいては量子コヒーレンスが必ずしも計算時間に亘って維持される必要がなく、基底状態が維持されれば良かったのはそのためである。アニールの概念はこのように幅広いものであり量子性の使い方に関しても量子アニールとは異なる方法論が有り得る。

その観点で発明されたのが前述の局所場応答法である。この方法では量子アニールと同様に、演算器としてのスピン系に時刻t = t₀で横磁場を印加し、磁場を徐々に縮小して時刻t = τで解を得る。演算器そのものは古典的であり、磁場に対してスピンが応答する際に量子力学的情報が付加される。この方法は古典的マシンで動作させるので室温動作が可能であり、極低温が必要な量子アニールの課題を解決する。特許文献１〜３や非特許文献１では経験的、あるいは類似の問題を解いた結果から量子効果を平均的に含む応答関数を決定し、量子効果を取り入れた。非特許文献２では量子力学の線形重ね合わせの性質を現象論的に取り入れることにより特許文献１〜３や非特許文献１の方法に比べて解精度を向上させた。しかし、量子力学のもうひとつの重要な性質である量子縺れの性質が十分に取り入れられたものではなかった。

WO2015/118639 WO2016/157333 WO2016/194221

T. Tomaru, "Quasi-Adiabatic Quantum Computing Treated with c-Numbers Using the Local-Field Response," J. Phys. Soc. Jpn. 85, 034802 (2016). T. Tomaru, "Improved Local Field Response Method Working as Quasi-Quantum Annealing," J. Phys. Soc. Jpn. 86, 054801 (2017).

以上述べたように、量子アニールは超伝導磁束量子ビットを用いるために極低温冷却装置を必要とする。また、局所場応答法は室温で動作するものの、量子力学の重要な性質である量子縺れの性質が十分に取り入れられたものではなく、性能が制限されていた。そこで本発明の目的は、全数探索を必要とするような難しい課題に対して十分な性能持った室温動作可能な計算機及び演算プログラムを提供することにある。

変数としてのスピンを局所的な有効磁場に応答させる局所場応答法において、時間軸を離散的とし、有効磁場に対してスピンを応答させる際に、スピンの大きさの小さいサイトから順番に有効磁場を求めてスピンを応答させ、スピンの符号が反転した場合はその情報を以降のサイトの有効磁場決定に反映させて、量子縺れに起因した多体効果を現象論的に取り入れる。より具体的には以下のようになる。

演算部、記憶部、制御部を具備し、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算機、あるいは計算機による計算方法であって、
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 ≦ τ)まで演算するものとする。ここでN、mは自然数である。

各時刻t_i(i ＝ 1, 2, .., m)で変数B_eff,j ^z(t_i)及びs_j ^z(t_i)をこの順番で定めるものとし、該B_eff,j ^z(t_i)はs_k ^z(t_i−1)、J_kj、g_j、t_iの関数であり、該s_j ^z(t_i)はB_eff,j ^z(t_i)及びt_iの関数であり、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、時刻t_i(i ＝ 1, 2, .., m)のB_eff,j ^z(t_i)及びs_j ^z(t_i)を求めるに当たっては、まずs_j ^z(t_i−1)を|s_m1 ^z(t_i−1)|≦|s_m2 ^z(t_i−1)|≦|s_m3 ^z(t_i−1)|≦・・・≦|s_mN ^z(t_i−1)|のように降べきの順に並べ、最初にサイトm₁のB_eff,m1 ^z(t_i)及びs_m1 ^z(t_i)を求めてs_m1 ^z(t_i−1) = sgn(s_m1 ^z(t_i))|s_m1 ^z(t_i−1)|とし、次にサイトm₂のB_eff,m2 ^z(t_i)及びs_m2 ^z(t_i)を求めてs_m2 ^z(t_i−1) = sgn(s_m2 ^z(t_i))|s_m2 ^z(t_i−1)|とし、続いてサイトm₃の演算を同様に実行し、以下サイトm_Nまで同様な演算を実行して時刻t_iの演算を行うものとし、
時刻ステップを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として解を定めることを特徴とする。

また、上記ではサイトm_Nまで同様な演算を実行して時刻t_iの演算を行っているが、m_x（ただし、0<x≦N）まで同様な演算を実行し、サイトm_x+1以降では、従来の局所場応答法と同様に、全てのサイトを独立かつ並列的に処理して時刻t_iの演算を行うように構成してもよい。なおxは自然数である。

本方法は古典的マシン上で動作させるものであり、極低温にする必要が無く、また量子コヒーレンスを考慮する必要もない。その結果、使用可能なリソースが広範囲になり電気回路等も利用できる。さらに、量子縺れの効果を現象論的に取り入れたことにより解精度が向上するとともに計算時間が短縮される。これらの性質により高い解精度で難問を解ける実用的な計算機が実現する。

実施例の原理を模式的に示した模式図である。 実施例１に係るアルゴリズムの一例をフローチャートで示した流れ図である。 応答関数r_bの具体的な値を示したグラフ図である。 実施例２に係るアルゴリズムの一例をフローチャートで示した流れ図である。 実施例３に係るアルゴリズムの一例をフローチャートで示した流れ図である。 実施例３に係るアルゴリズムの他の例をフローチャートで示した流れ図である。 実施例４に係るアルゴリズムの一例をフローチャートで示した流れ図である。 実施例５に係るアルゴリズムの一例をフローチャートで示した流れ図である。 実施例６に係るアルゴリズムの一例をフローチャートで示した流れ図である。 実施例７に係る計算機構成の一例を示したブロック図である。

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

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

同様なあるいは対応する機能を有する要素が複数ある場合には、同一の符号に異なる添字を付して説明する場合がある。ただし、複数の要素を区別する必要がない場合には、添字を省略して説明する場合がある。

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

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

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

(1)

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

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

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

(5)
のようにx, z成分だけで定義し、スピンの向きをこの応答関数に基づき定める。スピン系が古典的ならば各スピンの応答は各サイトの有効磁場だけで求まり、応答関数はr_b(t)=1になる。

しかし、量子力学には非局所相関(量子縺れ、entanglement)があり一般にr_b(t)≠1である。すでに言及したように式(5)は期待値をとることにより古典的な式に移行しているが、r_b(t)≠1を通して量子効果が取り込まれる。r_b(t)の値は経験的あるいは類似問題の量子力学的な事前計算により求める。経験的あるいは類似問題を元にするためにここでの量子効果は平均的なものになる。局所場応答法は様々な方法で量子効果を取り入れることが可能であるが、r_b(t)≠1を通した方法はその一つである。尚、局所場応答法は量子効果を取り入れずに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)は

(6)
となる。

局所場応答法の時間発展は離散的に行い、時刻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=t_mで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）。ここでt_m≦τとしたのはt=τまで時間発展させなくても解が収束することが多いからである。

手順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 _i )/B_eff, _j ^x (t _i ))}である。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に設定することは有効な手段である。

実施例１ではr_b(t)≠1を通して量子効果を平均的に取り入れられることを示した。しかし、量子効果は問題に依存すると共に各時刻で変化する。平均量だけでは十分に量子効果を取り込めない。本実施例では、量子縺れに係わる量子効果を各時刻のスピン状態に依存した形で現象論的に取り込む方法を示す。

量子縺れの影響は多体効果として現れる。量子縺れが大きければ，あるスピンが反転（符号が反転）する際に他のスピンが同時に反転する確率が大きくなる。図２のアルゴリズムではt = t_i−1のs_j ^z(t_i−1)を利用してt = t_iの有効磁場B_eff,j ^z(t_i)をサイトごとに求めた。サイトごとに独立した計算であり一体近似である。そのためスピンの同時反転が十分に考慮されていない。そこで図４のアルゴリズムによりスピンの同時反転を考慮する。

スピンが反転する近傍ではs_j ^z ≒ 0である。よって|s_j ^z|の値が小さいスピンほど反転する確率が高い。そこでまず各時刻t_i−1において|s_j ^z (t_i−1)|を降べきの順に並べる。

図４の111に示すように|s_j ^z (t_i−1)|が小さいサイトの順にm₁, m₂, m₃, ・・・とする。サイトm₁に対するB_eff,m1 ^z(t_i)を式(4)に従い求め（処理102a）、式(6)に従いs_m1 ^z (t_i)を求める（処理103）。s_m1 ^z(t_i−1)からs_m1 ^z (t_i)への時間発展で符号反転がなければサイトm₂のB_eff,m2 ^z(t_i)及びs_m2 ^z (t_i)の計算に進む。符号反転があればs_m1 ^z(t_i−1)を−s_m1 ^z(t_i−1)に変更してからサイトm₂のB_eff,m2 ^z(t_i)及びs_m2 ^z (t_i)の計算に進む。言い換えればs_m1 ^z(t_i−1) → sgn(s_m1 ^z(t_i))・|s_m1 ^z(t_i−1)|のようにs_m1 ^z(t_i−1)の符号をs_m1 ^z(t_i)の符号に付け替える（処理112）。なお、sgnは変数の符号に応じて1,−1のいずれかを返す符号関数を示す。

このようにすれば時刻t_iにおけるサイトm₁の変化がサイトm₂の時間発展に直ちに反映される。サイトm₁とm₂が量子縺れ状態になっていればサイトm₁がスピン反転すればサイトm₂もスピン反転する確率が高くなる。この効果がここで述べた符号の付け替えにより取り入れられる。即ち、量子縺れの効果が現象論的に取り入れた訳である。

s_m2 ^z(t_i)の計算結果もs_m1 ^z(t_i)と同様な符号処理を行う。即ち、s_m2 ^z(t_i−1) → sgn(s_m2 ^z(t_i))・|s_m2 ^z(t_i−1)|とする（処理112）。サイトm₃以降も同様な処理を繰り返してN個のs_j ^z(t_i)を求め（処理113）、時刻t_iの処理が完了する。

時刻t_iの処理が完了すれば時刻t_i+1の処理に進み、時刻t = t_m ≦τまで同様な処理を繰り返して最終的な解を得る。最終的な解は実施例１の場合と同様にs_j ^z(t_m)>0ならばs_j ^zfd=1、s_j ^z(t_m)<0ならばs_j ^zfd=−1である（処理201及び202）。

スピン反転は量子力学的にはトンネル現象が起こっているといった捉え方もできる。その捉え方に従えば、本実施例はMultiple tunnelingを考慮したものとも言える。

得られた解のエネルギー値は式(7)

(7)
により与えられる。局所場応答法は量子アニールと類似の動作をするものであり、t = 0で基底状態に準備したスピン系を時間発展させて理想的にはt =τで問題を設定した系の基底状態に導く。本実施例の方法は解の収束性が高く、t = t_mの状態は高い確率で演算中の最低エネルギー状態（高い確率で基底状態）である。しかし、実施例１の方法では解の収束性が悪く、最低エネルギー状態がt < t_mで現れることもあり、実施例１を用いた場合、各時刻で式(7)を用いてエネルギーを算出し、その中で最低エネルギーであった状態を最終解に選ぶ必要がある。そのための計算量は式(7)の第1項がO(N²)、第2項がO(N)である（Oはランダウの記号）。両者をまとめればO(N²)である。一方、本実施例の方法を用いれば解の収束性が高いのでエネルギーの計算をする必要がなくO(N²)の計算量を節約できる。

その一方で、本実施例の方法は処理111においてスピンの並べ替えを行う。この処理量は以下のように見積もれる。スピンごとに他のスピンと比較する処理を単純に繰り返せばO(N²)である。これが処理量の上限になる。下限は以下のように見積もれる。すでに降べきの順に並んでいるスピン列を考え，隣同士を比較して並べ替えるとする。並べ替えがなければ隣同士の比較だけなので計算量はO(N)である。一般に激しい並べ替えはめったにないので実際の計算量はO(N)に近くなる。よって本実施例のオーバーヘッドはO(N)程度になり実施例１のオーバーヘッドO(N²)に比べて小さい。

局所場応答法全体の計算量は処理102aにおける有効磁場の計算で決まり、Nサイトすべてに対してO(N)の計算をするのでO(N²)である。よってNが十分に大きい系では本実施例のオーバーヘッドO(N)程度は無視できることになる。

以上述べたように、図２の例では処理102a, 103が全てのサイトで独立かつ並列的に処理されていた。そのために処理速度は速いが量子縺れの考慮が不十分であった。一方、本実施例の図４の例では処理112, 111を加えることにより、同時刻の他のスピンの影響を処理102a, 103において反映させ、量子縺れの効果を現象論的に取り入れた。その結果、解精度が向上し、解の収束性が向上し、演算時間が節約できるようになった。

量子力学的には有効磁場は式(4)に基づき定まる。σ^_k ^zの固有値は±1である。しかし局所場応答法ではスピン変数s_k ^zが期待値<σ^_k ^z>の値を取るように動作させるので|s_k ^z|≦1である。従って、一般にg_jに比べてΣ_k（≠j）J_kjs_k ^zの項を過小評価することになる。

Σ_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)を局所項にすれば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となってしまうことに対処するためである。

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

図６には変形例を示している。パラメタc_aを導入し、因子c_iをc_a・c_iとすれば因子c_iの大きさを調整できる。その場合を図６に示す。後に説明する図９の処理10cも因子c_iをc_a・c_iとしている。c_aの値は経験的に決めるものであり1〜50程度である。

量子力学的スピン系では常にスピン同士で影響し合っている。即ち、あるサイトjのスピンσ^_j ^zは別のサイトkのスピンσ^_k ^zに影響し、逆にσ^_k ^zがσ^_j ^zに影響する。従って、スピンσ^_j ^zはサイトkのスピンσ^_k ^zを経由して自身に影響する。量子力学においてあるスピンの状態が相互作用相手のスピンの状態だけでなく自身のスピンの状態に依存するのはそのためである。相互作用を通した自身への影響の大きさはΣ_k（≠j）J_kj ²に比例する。実施例１では平均化された量子効果について述べた。Σ_k（≠j）J_kj ²は二乗の項なので平均しても残る。平均化した応答関数がr_b ⁰(t)≠1となったのはこのためである。図３の円形の記号及び実線はその様子を示す。Σ_k（≠j）J_kj ²に関する詳細な理論は非特許文献２に記されている。

r_b ⁰(t)≠1の原因となったΣ_k（≠j）J_kj ²には問題ごとの情報J_kjが含まれている。この情報を利用できれば演算がより正確になる。そこで平均的応答関数r_b ⁰(t)をr_b ⁰ _mod(t)に置き換える。r_b ⁰ _mod(t)は1−r_b ⁰ _mod(t)=(1−r_b ⁰(t))Σ_k（≠j）J_kj ²/Σ_k（≠j）ave(J_kj ²)で定義される。ここでave(J_kj ²)はr_b ⁰(t)を決定する際に使用した問題のJ_kj ²を平均したものである。この改良により応答関数が実際の問題を反映したものになり解精度が向上する。

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

量子力学の特徴的な性質として量子縺れと線形重ね合わせが挙げられる。前者の量子縺れの効果は実施例４までに現象論的に取り入れられた。後者の線形重ね合わせの効果は非特許文献２において現象論的に取り入れられている。本実施例では両者を同時に取り入れる。その場合を示したのが図８である。

図８に示す実施例において、図７との違いは手順102bが手順102dに変わることである。

線形重ね合わせが顕著になるのはスピンの符号が変わる時間帯である。その近傍では量子力学的にはバンドが反交差となり、線形重ね合わせ状態となってs_j ^z(t)≒0であり、またそれに連動してB_eff,j ^z(t)≒0になる。この効果を現象論的に取り入れるためには時刻t_iにおける有効磁場B_eff,j ^z(t_i)を、t=t_iとそのひとつ前のt=t_i−1における有効磁場の線形結合にすればよい。具体的には時刻t_i−1におけるスピンの値s_j ^z(t_i−1)により式(8)を使って

(8)
B_j ^z0(t_i)を求める。次にひとつ前の時刻t_i−1における因子も考慮して式(9)により

(9)
B_j ^z(t_i)を求める。ここでuは0≦u≦1で解精度が高くなるように適当に定める。典型的な値はu≒0.1である。横磁場及びそのスケジュールも含めて有効磁場を記述すれば式(10)になる。

(10)
有効磁場が求まれば応答関数r_b ⁰ _mod(t)を利用して

(11)
式(11)からs_j ^z(t)を求める。

実施例５まで、量子縺れの効果を中心に各種量子効果を加えて解精度を向上させる実施例を述べた。しかし、以上の方法を用いても必ず正解に辿り着くとは限らない。そこで本実施例では補助的な手段を述べる。

局所場応答法は決定論的な方法なので同じスピンの初期値で同じ処理をすれば結果はいつも同じである。初期値や処理過程を変更すれば結果が異なる場合もある。そこで初期値や処理過程を変えて複数回磁場挿引の処理を実施し、その中で最低エネルギーを与えた場合のスピン状態を最終解に選んで解精度をさらに向上させる。その場合のフローチャートが図９である。

図９において、処理10a - 10dは図２や図４〜８で示した処理である。処理10a - 10dで得られたスピン値s_j ^zfdからそれぞれエネルギーH_p ^q(t_m) (q = 1, 2, ・・・)を処理303で求め、その最低値E_bestを処理304で求める。最低エネルギーを与えたスピン値が最終的な解になる。

処理10aは図２や図４〜８で示した処理をそのまま実行する。処理10bでは処理10aで得られた結果s_j ^zfdの符号反転したものを初期値とする。正確にはパラメタdivを用いて初期値をs_j ^z(t₀) = s_j ^zfd/divにする。divで割るのは初期値を十分に小さくするためで、divはm程度の大きさにする。符号反転するのは符号反転した状態がスピンの配位空間上、反転前の状態から最も遠い点であり、局所解に陥った場合に抜け出す機会を与えることができるからである。処理10cでは10aと同様にスピンの初期値を0にするが、局所項g_jに掛かる係数c_i（実施例３参照）に因子c_aを掛けてc_a・c_iとすることにより相互作用項と局所項のバランスを変更する。c_aの値は経験的に決めるものであり1〜50程度である。処理10dはスピンの初期値を乱数で決める。乱数を使うことにより多くの可能性を探れると共に処理を単純化できる。処理10dは乱数の初期値を変えて繰り返し実行する。実行回数が多いほど解精度は向上する。

本実施例はアルゴリズムとして示されており、通常のコンピュータ上でソフトウェアとして動作させることも専用のハードウェア上で動作させることもできる。本実施例の局所場応答法は比較的単純な演算を繰り返すのが特徴である。従って、繰り返し演算の部分を専用ハードウェアで構築し、その他の部分を汎用装置で実現するのが効果的である。

図１０に本実施例の計算機構成の一例を示す。図１０は通常の計算機の構成と類似であるが局所場応答演算装置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である。

実施例６では類似の処理10a - 10dを実施し、処理ごとに解s_j ^zfdを得た。得られた解は局所場応答演算装置600から主記憶装置501にデータ転送し、一般演算装置502を利用してエネルギーH_p ^q(t_m)やE_bestの計算を行う。即ち、繰り返し演算に属さない個別の演算は一般演算装置502を利用して局所場応答演算装置600の専用性を高める。

組合せ最適化問題の中で困難度の高い問題は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>jw_ijs_is_jを最小化する問題となる。イジング・スピングラスのハミルトニアンは

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

図４に示した実施例２では、スピンの数であるNサイトの変数全てについて量子縺れを考慮した計算を行っている。しかし、|s_j ^z|の値が大きくなるとスピンが反転する確率が低くなるので、所定以降のスピンについては変化がないものとして量子縺れの効果を無視してもよい。すなわち、図４の処理113をl=x<Nの条件に変更して、Yesの場合にはx+1番目以降のスピンについては処理112を行わずに、処理102a, 103を全てのサイトで独立かつ並列的に処理して時刻t_iにおける演算をするように変更しても良い。このようにすると、ある程度の精度を確保しつつ、処理時間を短縮することが可能である。

本実施例中、ソフトウェアで構成した機能と同等の機能は、FPGA（Field Programmable Gate Array）、ASIC（Application Specific Integrated Circuit）などのハードウェアでも実現できる。

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

10 - 304 フローチャートにおける各手順を表す
501 主記憶装置
502 一般演算装置
503 制御装置
504 補助記憶装置
505 入力装置
506 出力装置
600 局所場応答演算装置

Claims (13)

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，．．， m)で変数B_eff，j ^z(t_i)及びs_j ^z(t_i)をこの順番で定めるものとし、該B_eff，j ^z(t_i)はs_k ^z(t_i−1)、J_kj、g_j、t_iの関数であり、該s_j ^z(t_i)はB_eff，j ^z(t_i)及びt_iの関数であり、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、時刻t_i(i ＝ 1，2，．．，m)のB_eff，j ^z(t_i)及びs_j ^z(t_i)を求めるに当たっては、まずs_j ^z(t_i−1)を|s_m1 ^z(t_i−1)|≦|s_m2 ^z(t_i−1)|≦|s_m3 ^z(t_i−1)|≦…≦|s_mN ^z(t_i−1)|のように降べきの順に並べ、最初にサイトm₁のB_eff，m1 ^z(t_i)及びs_m1 ^z(t_i)を求めてs_m1 ^z(t_i−1) = sgn(s_m1 ^z(t_i))|s_m1 ^z(t_i−1)|とし、次にサイトm₂のB_eff，m2 ^z(t_i)及びs_m2 ^z(t_i)を求めてs_m2 ^z(t_i−1) = sgn(s_m2 ^z(t_i))|s_m2 ^z(t_i−1)|とし、続いてサイトm₃の演算を同様に実行し、以下サイトm_x(ただし3≦x≦N)まで同様な演算を実行して時刻t_iの演算を行うものとし、
   時刻ステップを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として解を定めることを特徴とする計算機。
2. 前記B_eff，j ^z(t_i)はB_eff，j ^z(t_i) = (t_i/τ)・(Σ_k（≠j）J_kjs_k ^z(t_i−1) ＋ g_j)により定めることを特徴とする請求項1記載の計算機。
3. ある定数γを用いて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を使ってs_j ^z(t_i) = f(B_eff，j ^z(t_i)，t_i) = sin{arctan(B_eff，j ^z(t_i)/B_eff，j ^x(t_i))}となることを特徴とする請求項１記載の計算機。
4. 前記関数fに関して補正パラメタr_s (t)及びr_b (t)を追加し、
   tanθ = r_b (t)・B_eff，j ^z(t_i)/B_eff，j ^x(t_i)によりθを定義し、s_j ^z(t_i) = r_s (t)・sinθによって前記s_j ^z(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_i)/B_eff，j ^x(t_i))}となることを特徴とする請求項３記載の計算機。
5. c_i = (Σ_k(s_k ^z(t_i−1))²/N)^1/2とし、g_j ^norm(t_i) = c_i・g_jとして、前記B_eff，j ^z(t_i)をB_eff，j ^z(t_i) = (t_i/τ)・(Σ_k（≠j）J_kjs_k ^z (t_i−1)＋ g_j ^norm(t_i))により定めることを特徴とする請求項２記載の計算機。
6. あるパラメタc_aを利用して前記B_eff，j ^z(t_i)をB_eff，j ^z(t_i) = (t_i/τ)・(Σ_k（≠j）J_kjs_k ^z (t_i−1)＋ c_a・g_j ^norm(t_i))により定めることを特徴とする請求項５記載の計算機。
7. 前記補正パラメタr_b (t)に対してδr_b (t)≡1−r_b (t)を定義し、δr_b(t)∝Σ_k（≠j）J_kj ²とすることを特徴とする請求項４記載の計算機。
8. B_j ^z0(t_i) = (Σ_k（≠j）J_kjs_k ^z(t_i−1)＋ g_j ^norm(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_eff，j ^z(t_i)をB_eff，j ^z(t_i) = B_j ^z(t_i)・t_i/τにより定めることを特徴とする請求項５記載の計算機。
9. s_j ^zfdを求める請求項1記載の演算を複数回実行するものとし、パラメタdivを前記m程度の大きさの値とし、2回目以降の演算の初期値を1回目以前の解s_j ^zfdを利用してs_j ^z(t₀) =−s_j ^zfd/divとするか、乱数を利用してs_j ^z(t₀) = 1/divあるいはs_j ^z(t₀) = −1/divとし、それぞれの演算ごとにH_p ＝−Σ_k＞jJ_kjs_k ^zfd(t_i)s_j ^zd −Σ_jg_js_j ^zfdを算出し、H_pが最小値となった演算のs_j ^zfdを最終解とすることを特徴とする請求項１記載の計算機。
10. サイトm_x＋1以降は、全てのサイトを独立かつ並列的に処理して時刻t_iの演算を行う、
    請求項１記載の計算機。
11. 演算部、記憶部、制御部を具備する計算機を用い、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算方法であって、
    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，．．，m)で変数B_eff，j ^z(t_i)及びs_j ^z(t_i)をこの順番で定めるものとし、
    該B_eff，j ^z(t_i)はs_k ^z(t_i−1)、J_kj、g_j、t_iの関数であり、該s_j ^z(t_i)はB_eff，j ^z(t_i)及びt_iの関数であり、時刻t₀の初期値はB_j ^z(t₀)=0及びs_j ^z(t₀)=0とし、
    時刻t_i(i ＝ 1， 2， ．．， m)のB_eff，j ^z(t_i)及びs_j ^z(t_i)を求めるに当たっては、まずs_j ^z(t_i−1)を|s_m1 ^z(t_i−1)|≦|s_m2 ^z(t_i−1)|≦|s_m3 ^z(t_i−1)|≦…≦|s_mN ^z(t_i−1)|のように降べきの順に並べ、
    最初にサイトm₁のB_eff，m1 ^z(t_i)及びs_m1 ^z(t_i)を求めてs_m1 ^z(t_i−1) = sgn(s_m1 ^z(t_i))|s_m1 ^z(t_i−1)|とし、それ以降も同様の演算をサイトm_x(ただし1≦x≦N)まで実行して時刻t_iの演算を行うものとし、
    時刻ステップを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として解を定めることを特徴とする計算方法。
12. サイトm₂以降も同様の演算をサイトm_Nまで実行して時刻t_iの演算を行う、請求項１１記載の計算方法。
13. サイトm_x＋1以降は、全てのサイトを独立かつ並列的に処理して時刻t_iの演算を行う、請求項１１記載の計算方法。

Images