(57)【要約】【課題】全数探索を必要とするような組み合わせ最適化問題に対して十分な性能持った室温動作可能な計算技術を提供することである。【解決手段】  変数としてのスピンを局所的な有効磁場に応答させる局所場応答法において、時間軸を離散的とし、有効磁場に対してスピンを応答させる際に、スピンの大きさの小さいサイトから順番に有効磁場を求めてスピンを応答させ、スピンの符号が反転した場合はその情報を以降のサイトの有効磁場決定に反映させて、量子縺れに起因した多体効果を現象論的に取り入れる。【選択図】  図４
2019079225.tif 000002
【特許請求の範囲】
【請求項１】
演算部、記憶部、制御部を具備し、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算機であって、
N個の変数sjz (j = 1， 2，…，N)が−1≦sjz≦1の値域を取り、局所項を表す係数gjと変数間相互作用を表す係数Jkj (k， j ＝ 1， 2，…，N)によって課題の設定を行い、
前記演算部では、時刻をm分割して離散的にt = t0 (t0 = 0)からtm (tm ≦ τ)まで演算するものとし、
各時刻ti(i ＝ 1，2，．．， m)で変数Beff，jz(ti)及びsjz(ti)をこの順番で定めるものとし、該Beff，jz(ti)はskz(ti−1)、Jkj、gj、tiの関数であり、該sjz(ti)はBeff，jz(ti)及びtiの関数であり、時刻t0の初期値はBjz(t0)=0及びsjz(t0)=0とし、時刻ti(i ＝ 1，2，．．，m)のBeff，jz(ti)及びsjz(ti)を求めるに当たっては、まずsjz(ti−1)を|sm1z(ti−1)|≦|sm2z(ti−1)|≦|sm3z(ti−1)|≦…≦|smNz(ti−1)|のように降べきの順に並べ、最初にサイトm1のBeff，m1z(ti)及びsm1z(ti)を求めてsm1z(ti−1) = sgn(sm1z(ti))|sm1z(ti−1)|とし、次にサイトm2のBeff，m2z(ti)及びsm2z(ti)を求めてsm2z(ti−1) = sgn(sm2z(ti))|sm2z(ti−1)|とし、続いてサイトm3の演算を同様に実行し、以下サイトmx(ただし3≦x≦N)まで同様な演算を実行して時刻tiの演算を行うものとし、
時刻ステップをt = t0からt = tmに進めるにつれて前記変数sjzを−1あるいは1に近づけ、最終的にsjz < 0ならばsjzfd =−1、sjz ＞0ならばsjzfd = 1として解を定めることを特徴とする計算機。
【請求項２】
前記Beff，jz(ti)はBeff，jz(ti) = (ti/τ)・(Σk（≠j）Jkjskz(ti−1) ＋ gj)により定めることを特徴とする請求項1記載の計算機。
【請求項３】
ある定数γを用いてBeff，jx(ti) = γ(1 − ti/τ)としてtanθ = Beff，jz(ti)/Beff，jx(ti)によりθを定義し、前記sjz(ti)をsjz(ti) = sinθによって定めることとし、従って関数fを使ってsjz(ti) = f(Beff，jz(ti)，ti) = sin{arctan(Beff，jz(ti)/Beff，jx(ti))}となることを特徴とする請求項１記載の計算機。
【請求項４】
前記関数fに関して補正パラメタrs及びrbを追加し、
tanθ = rb・Beff，jz(ti)/Beff，jx(ti)によりθを定義し、sjz(ti) = rs・sinθによって前記sjz(ti)を定めることとし、従って前記関数fがf(Beff，jz(ti)， ti) ＝ rs・sin{arctan(rb・Beff，jz(ti)/Beff，jx(ti))}となることを特徴とする請求項３記載の計算機。
【請求項５】
ci = (Σk(skz(ti−1))2/N)1/2とし、gjnorm(ti) = ci・gjとして、前記Beff，jz(ti)をBeff，jz(ti) = (ti/τ)・(Σk（≠j）Jkjsk(ti−1)＋ gjnorm(ti))により定めることを特徴とする請求項２記載の計算機。
【請求項６】
あるパラメタcaを利用して前記Beff，jz(ti)をBeff，jz(ti) = (ti/τ)・(Σk（≠j）Jkjsk(ti−1)＋ ca・gjnorm(ti))により定めることを特徴とする請求項５記載の計算機。
【請求項７】
前記補正パラメタrbに対してδrb≡1−rbを定義し、δrb(t)∝Σk（≠j）Jkj2とすることを特徴とする請求項４記載の計算機。
【請求項８】
Bjz0(ti) = (Σk（≠j）Jkjskz(ti−1)＋ gjnorm(ti))を定義し、0≦u≦1を満たすパラメタuを用いてBjz(ti) = (1−u)Bjz0(ti)＋ uBjz(ti−1)を定義し、前記Beff，jz(ti)をBeff，jz(ti) = Bjz(ti)・ti/τにより定めることを特徴とする請求項５記載の計算機。
【請求項９】
sjzfdを求める請求項1記載の演算を複数回実行するものとし、パラメタdivを前記m程度の大きさの値とし、2回目以降の演算の初期値を1回目以前の解sjzfdを利用してsjz(t0) =−sjzfd/divとするか、乱数を利用してsjz(t0) = 1/divあるいはsjz(t0) = −1/divとし、それぞれの演算ごとにHp ＝−Σk＞jJkjskzfd(ti)sjzd −Σjgjsjzfdを算出し、Hpが最小値となった演算のsjzfdを最終解とすることを特徴とする請求項１記載の計算機。
【請求項１０】
サイトmx＋1以降は、全てのサイトを独立かつ並列的に処理して時刻tiの演算を行う、請求項１記載の計算機。
【請求項１１】
演算部、記憶部、制御部を具備する計算機を用い、前記制御部の制御により、前記記憶部と前記演算部との間でデータをやり取りしながら演算を行う計算方法であって、
N個の変数sjz (j = 1，2，…， N)が−1≦sjz≦1の値域を取り、局所項を表す係数gjと変数間相互作用を表す係数Jkj (k， j ＝ 1， 2，…，N)によって課題の設定を行い、
前記演算部では、時刻をm分割して離散的にt = t0 (t0 = 0)からtm (tm ≦ τ)まで演算するものとし、
各時刻ti(i ＝ 1，2，．．，m)で変数Beff，jz(ti)及びsjz(ti)をこの順番で定めるものとし、
該Beff，jz(ti)はskz(ti−1)、Jkj、gj、tiの関数であり、該sjz(ti)はBeff，jz(ti)及びtiの関数であり、時刻t0の初期値はBjz(t0)=0及びsjz(t0)=0とし、
時刻ti(i ＝ 1， 2， ．．， m)のBeff，jz(ti)及びsjz(ti)を求めるに当たっては、まずsjz(ti−1)を|sm1z(ti−1)|≦|sm2z(ti−1)|≦|sm3z(ti−1)|≦…≦|smNz(ti−1)|のように降べきの順に並べ、
最初にサイトm1のBeff，m1z(ti)及びsm1z(ti)を求めてsm1z(ti−1) = sgn(sm1z(ti))|sm1z(ti−1)|とし、それ以降も同様の演算をサイトmx(ただし1≦x≦N)まで実行して時刻tiの演算を行うものとし、
時刻ステップをt = t0からt = tmに進めるにつれて前記変数sjzを−1あるいは1に近づけ、最終的にsjz < 0ならばsjzfd = −1、sjz ＞ 0ならばsjzfd = 1として解を定めることを特徴とする計算方法。
【請求項１２】
サイトm2以降も同様の演算をサイトmNまで実行して時刻tiの演算を行う、請求項１１記載の計算方法。
【請求項１３】
サイトmx＋1以降は、全てのサイトを独立かつ並列的に処理して時刻tiの演算を行う、請求項１１記載の計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、全数探索を必要とするような組み合わせ最適化問題等に対して高速演算を可能にする計算技術に関するものである。
【背景技術】
【０００２】
IoT (Internet of Things)といった言葉に代表されるように、現代はあらゆるモノがインターネットにつながり，モノから情報が収集され，総合した情報を元にモノが制御される。制御では多くの選択肢の中から最適解が見出されて実行される。極論すれば現代の情報技術は最適解を見出すものとも言える。
【０００３】
こういった背景の中で近年注目されるようになってきたのが量子アニール、別名断熱量子計算と呼ばれる手法である。この方法は、ある物理系の基底状態が解になるように問題を設定し基底状態を見つけることを通して解を得ようとするものである。問題を設定した物理系のハミルトニアンをH^pとする。但し、演算開始時のハミルトニアンはH^pではなく、基底状態に準備しやすい別のハミルトニアンH^0とする。次に十分に時間を掛けてハミルトニアンをH^0からH^pに移行させる。十分に時間を掛ければ系は基底状態に居続け、最終的にハミルトニアンH^pの基底状態（解状態）を得る。これが量子アニールの原理である。
【０００４】
イジング・スピングラスと呼ばれる物理系を利用した基底状態探索法はNP困難と呼ばれる問題にも対応できる。また組合せ最適化問題の中で困難度の高い問題はNP困難に属する。さらに計算複雑性理論でPに分類される問題やNPに分類される問題はすべてNP困難問題に帰着できる。よって、イジング・スピングラス系で量子アニールを適用すれば組合せ最適化問題をほぼすべて解けることになり、情報技術の最も重要な課題が解決される。
【０００５】
量子アニールが注目されるもう一つの理由はディコヒーレンスに対して頑強なことである。量子コンピュータでは量子コヒーレンスが計算時間に亘って保たれていなければならなかったが、量子アニールでは条件が緩和され基底状態が維持されていれば正解が得られる。必ずしも量子コヒーレンスが維持されている必要はない。現状の技術レベルでは純粋な量子系を構築することは困難であり、よって量子コヒーレンスを計算時間に亘って維持することは困難である。この点を考慮すれば量子アニールが注目される理由が理解できる。但し、量子アニールにも欠点がある。量子アニールを実現しうるのは現状では超伝導磁束量子ビット系に限られており、極低温冷却装置を必要とする。極低温の必要性は実用的なコンピュータ実現のためには課題である。
【０００６】
この課題を解決するために考案された方法が以下で述べる局所場応答法である（特許文献１〜３、非特許文献１、２）。まず量子アニールを再考する。アニール（焼きなまし）の概念は元々量子・古典に関係なく存在するものであり、量子アニールは量子性を使って古典アニールの性能を向上させようとしたものである。量子アニールにおいては量子コヒーレンスが必ずしも計算時間に亘って維持される必要がなく、基底状態が維持されれば良かったのはそのためである。アニールの概念はこのように幅広いものであり量子性の使い方に関しても量子アニールとは異なる方法論が有り得る。
【０００７】
その観点で発明されたのが前述の局所場応答法である。この方法では量子アニールと同様に、演算器としてのスピン系に時刻t = t0で横磁場を印加し、磁場を徐々に縮小して時刻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個の変数sjz (j = 1, 2,・・・ , N)が−1≦sjz≦1の値域を取り、局所項を表す係数gjと変数間相互作用を表す係数Jkj (k, j ＝ 1, 2, ・・・, N)によって課題の設定を行い、
前記演算部では、時刻をm分割して離散的にt = t0 (t0 = 0)からtm (tm ≦ τ)まで演算するものとする。ここでN、mは自然数である。
【００１３】
各時刻ti(i ＝ 1, 2, .., m)で変数Beff,jz(ti)及びsjz(ti)をこの順番で定めるものとし、該Beff,jz(ti)はskz(ti−1)、Jkj、gj、tiの関数であり、該sjz(ti)はBeff,jz(ti)及びtiの関数であり、時刻t0の初期値はBjz(t0)=0及びsjz(t0)=0とし、時刻ti(i ＝ 1, 2, .., m)のBeff,jz(ti)及びsjz(ti)を求めるに当たっては、まずsjz(ti−1)を|sm1z(ti−1)|≦|sm2z(ti−1)|≦|sm3z(ti−1)|≦・・・≦|smNz(ti−1)|のように降べきの順に並べ、最初にサイトm1のBeff,m1z(ti)及びsm1z(ti)を求めてsm1z(ti−1) = sgn(sm1z(ti))|sm1z(ti−1)|とし、次にサイトm2のBeff,m2z(ti)及びsm2z(ti)を求めてsm2z(ti−1) = sgn(sm2z(ti))|sm2z(ti−1)|とし、続いてサイトm3の演算を同様に実行し、以下サイトmNまで同様な演算を実行して時刻tiの演算を行うものとし、
時刻ステップをt = t0からt = tmに進めるにつれて前記変数sjzを−1あるいは1に近づけ、最終的にsjz < 0ならばsjzfd = −1、sjz > 0ならばsjzfd = 1として解を定めることを特徴とする。
【００１４】
また、上記ではサイトmNまで同様な演算を実行して時刻tiの演算を行っているが、mx（ただし、0<x≦N）まで同様な演算を実行し、サイトmx+1以降では、従来の局所場応答法と同様に、全てのサイトを独立かつ並列的に処理して時刻tiの演算を行うように構成してもよい。なおxは自然数である。
【発明の効果】
【００１５】
本方法は古典的マシン上で動作させるものであり、極低温にする必要が無く、また量子コヒーレンスを考慮する必要もない。その結果、使用可能なリソースが広範囲になり電気回路等も利用できる。さらに、量子縺れの効果を現象論的に取り入れたことにより解精度が向上するとともに計算時間が短縮される。これらの性質により高い解精度で難問を解ける実用的な計算機が実現する。
【図面の簡単な説明】
【００１６】
【図１】実施例の原理を模式的に示した模式図である。
【図２】実施例１に係るアルゴリズムの一例をフローチャートで示した流れ図である。
【図３】応答関数rbの具体的な値を示したグラフ図である。
【図４】実施例２に係るアルゴリズムの一例をフローチャートで示した流れ図である。
【図５】実施例３に係るアルゴリズムの一例をフローチャートで示した流れ図である。
【図６】実施例３に係るアルゴリズムの他の例をフローチャートで示した流れ図である。
【図７】実施例４に係るアルゴリズムの一例をフローチャートで示した流れ図である。
【図８】実施例５に係るアルゴリズムの一例をフローチャートで示した流れ図である。
【図９】実施例６に係るアルゴリズムの一例をフローチャートで示した流れ図である。
【図１０】実施例７に係る計算機構成の一例を示したブロック図である。
【発明を実施するための形態】
【００１７】
実施の形態について、図面を用いて詳細に説明する。ただし、本発明は以下に示す実施の形態の記載内容に限定して解釈されるものではない。本発明の思想ないし趣旨から逸脱しない範囲で、その具体的構成を変更し得ることは当業者であれば容易に理解される。
【００１８】
以下に説明する発明の構成において、同一部分又は同様な機能を有する部分には同一の符号を異なる図面間で共通して用い、重複する説明は省略することがある。
【００１９】
同様なあるいは対応する機能を有する要素が複数ある場合には、同一の符号に異なる添字を付して説明する場合がある。ただし、複数の要素を区別する必要がない場合には、添字を省略して説明する場合がある。
【実施例１】
【００２０】
実施例１では量子力学的な記述から出発し、それを古典的な形式に移行することを通して土台となる原理を述べる。
【００２１】
図１に本実施例の原理を模式的に示す。基本的枠組みは特許文献１及び非特許文献１に記載の局所場応答法と同じである。t = 0において横磁場を印加してスピンを一方向に揃える。その後、横磁場をゆっくりと減少させてt = τで問題設定のハミルトニアンにする。スピンは各時刻でそれぞれに掛かる局所的有効磁場に応答して時間発展する。
【００２２】
問題設定のハミルトニアンとt = 0におけるハミルトニアンをそれぞれ
【００２３】
【数１】
2019079225.tif 000003
【００２４】
(1)
【００２５】
【数２】
2019079225.tif 000004
【００２６】
(2)
とし、時刻tにおけるハミルトニアンを
【００２７】
【数３】
2019079225.tif 000005
【００２８】
(3)
とする。τが演算時間である。１スピン系の類推からサイトjのスピンが受ける有効磁場はB^eff,j = −∂H^/∂σ^jで与えられる。
【００２９】
【数４】
2019079225.tif 000006
【００３０】
(4)
本実施例の局所場応答法は期待値を取った<σ^j>をスピン変数とみなして古典的マシン上で動作させるものである。式(1)、 (2)から明らかなように<σ ^j>及び<B^eff,j>はx, z成分だけからなる。そこで応答関数rb(t)を
【００３１】
【数５】
2019079225.tif 000007
【００３２】
(5)
のようにx, z成分だけで定義し、スピンの向きをこの応答関数に基づき定める。スピン系が古典的ならば各スピンの応答は各サイトの有効磁場だけで求まり、応答関数はrb(t)=1になる。
【００３３】
しかし、量子力学には非局所相関(量子縺れ、entanglement)があり一般にrb(t)≠1である。すでに言及したように式(5)は期待値をとることにより古典的な式に移行しているが、rb(t)≠1を通して量子効果が取り込まれる。rb(t)の値は経験的あるいは類似問題の量子力学的な事前計算により求める。経験的あるいは類似問題を元にするためにここでの量子効果は平均的なものになる。局所場応答法は様々な方法で量子効果を取り入れることが可能であるが、rb(t)≠1を通した方法はその一つである。尚、局所場応答法は量子効果を取り入れずにrb(t)=1として動作させても良い。量子効果を含めなくても局所場応答法自体は動作する。
【００３４】
式(5)の４つの変数<σ^jz(ti)>、<σ^jx(ti)>、<B^eff,jz(ti)>、<B^eff,jx(ti)>は期待値を取っており古典的な量である。そこで量子力学の記法から古典物理の記法に変更する。即ち、<σ^jx(ti)>→sjx(ti) <σ^jz(ti)>→sjz(ti)、<B^eff,jx(ti)>→Beff,jx(ti) <B^eff,jz(ti)>→Beff,jz(ti)とする。この記法の変更により式(5)は
【００３５】
【数６】
2019079225.tif 000008
【００３６】
(6)
となる。
【００３７】
局所場応答法の時間発展は離散的に行い、時刻tiにおけるBeff,jz(ti)は式(4)に従い時刻ti−1におけるskz(ti−1)から決める。時刻tiにおけるsjz(ti)は式(6)に従い時刻tiにおけるBeff,jz(ti)を元に決める。この手続きを繰り返す。図２はこれをフローチャートとしてまとめたものである。
【００３８】
図２の101はアルゴリズムの出発点で初期値の設定を表す。102aでは時刻ti−1におけるskz(ti−1)を使ってBeff,jz(ti)を求める。また時刻に依存して横磁場強度Beff,jx(ti)を定める。103ではBeff,jz(ti)/Beff,jx(ti)と応答関数rb(t)によりスピンの向きに相当するtanθ= rb(t)・Beff,jz(ti)/Beff,jx(ti)を求め、スピンの大きさを表すパラメタrs(t)を使ってsjz(ti) = rs(t)・sinθを求める。rs(t)はrs(t)2= sjx(ti)2+ sjz(ti)2で定義される量で、rb(t)と同様に経験的あるいは事前計算で決めておく。具体例は２つ後のパラグラフで述べる。図２の102aと103が繰り返し計算の１セットである。このセットをt=tm(≦τ)まで繰り返し、t=tmでsjz(tm)>0ならばsjzfd=1、sjz(tm)<0ならばsjzfd=−1として解sjzfdを得る（処理201及び202）。ここでtm≦τとしたのはt=τまで時間発展させなくても解が収束することが多いからである。
【００３９】
手順103はsjz(ti) = f(Beff,jz(ti),ti)のように関数fを使って一般的に書くこともでき、f(Beff,jz(ti),ti) = rs(t)・sin{arctan(rb(t)・Beff,jz(tk)/Beff,jx(tk))}である。rs(t)はスピンの大きさを表すパラメタで0≦rs(t)≦1である。また、−1≦sjz(ti)≦1である。ここで、rb(t)=1及びrs(t)=1ならば純粋に古典的になる。
【００４０】
図３に応答関数rb(t)の事前計算の一例を示す。これは８ビット系でJij及びgjを[−5、5]の一様乱数で決めた場合である。100個の問題の結果であり800点からなる（100個の問題×８ビット）。ここで、Bzx≡<B^eff,jz>/<B^eff,jx>、szx≡<σ^jz>/<σ^jx>である。点が厳密に量子力学的に求めた値である。量子力学の非局所相関を反映して応答関数が大きくばらつく。丸は横軸を４０分割して平均を取ったものである。平均化されて応答関数は滑らかなBzx依存性になる。滑らかならば数個のパラメタで記述可能になる。非特許文献１では４つのパラメタを使って滑らかな応答関数を記述する手法を述べており、図３の実線rb0(t)はその手法により求めたものである。もうひとつのパラメタrs(t)も同じ４つのパラメタをから求まる。
【００４１】
以上、図３において応答関数の一例を示し、非特許文献１を引用してrb0(t)及びrs(t)の決定法について言及した。しかし、応答関数rb0(t)やrs(t)の決定法はそれに制限されるものではなく様々な方法がありえると共に経験的に定めることもできる。以下の実施例ではrb0(t)を図３の実線に限定せずに量子効果を平均的に取り入れた応答関数を表すものとする。
【００４２】
尚、rb(t)と同様にrs(t)もrs(t)=1として動作させることも可能である。rb(t)の値域が−∞<rb(t)<∞であるのに対して0≦rs(t)≦1なので、rs(t)=1としたことによる最終解に対する影響はrb(t)=1としたことによる影響に比べて小さい。従って、rs(t)決定のための事前情報が不十分な場合にrs(t)=1に設定することは有効な手段である。
【実施例２】
【００４３】
実施例１ではrb(t)≠1を通して量子効果を平均的に取り入れられることを示した。しかし、量子効果は問題に依存すると共に各時刻で変化する。平均量だけでは十分に量子効果を取り込めない。本実施例では、量子縺れに係わる量子効果を各時刻のスピン状態に依存した形で現象論的に取り込む方法を示す。
【００４４】
量子縺れの影響は多体効果として現れる。量子縺れが大きければ，あるスピンが反転（符号が反転）する際に他のスピンが同時に反転する確率が大きくなる。図２のアルゴリズムではt = ti−1のsjz(ti−1)を利用してt = tiの有効磁場Beff,jz(ti)をサイトごとに求めた。サイトごとに独立した計算であり一体近似である。そのためスピンの同時反転が十分に考慮されていない。そこで図４のアルゴリズムによりスピンの同時反転を考慮する。
【００４５】
スピンが反転する近傍ではsjz ≒ 0である。よって|sjz|の値が小さいスピンほど反転する確率が高い。そこでまず各時刻ti−1において|sjz (ti−1)|を降べきの順に並べる。
【００４６】
図４の111に示すように|sjz (ti−1)|が小さいサイトの順にm1, m2, m3, ・・・とする。サイトm1に対するBeff,m1z(ti)を式(4)に従い求め（処理102a）、式(6)に従いsm1z (ti)を求める（処理103）。sm1z(ti−1)からsm1z (ti)への時間発展で符号反転がなければサイトm2のBeff,m2z(ti)及びsm2z (ti)の計算に進む。符号反転があればsm1z(ti−1)を−sm1z(ti−1)に変更してからサイトm2のBeff,m2z(ti)及びsm2z (ti)の計算に進む。言い換えればsm1z(ti−1) → sgn(sm1z(ti))・|sm1z(ti−1)|のようにsm1z(ti−1)の符号をsm1z(ti)の符号に付け替える（処理112）。なお、sgnは変数の符号に応じて1,−1のいずれかを返す符号関数を示す。
【００４７】
このようにすれば時刻tiにおけるサイトm1の変化がサイトm2の時間発展に直ちに反映される。サイトm1とm2が量子縺れ状態になっていればサイトm1がスピン反転すればサイトm2もスピン反転する確率が高くなる。この効果がここで述べた符号の付け替えにより取り入れられる。即ち、量子縺れの効果が現象論的に取り入れた訳である。
【００４８】
sm2z(ti)の計算結果もsm1z(ti)と同様な符号処理を行う。即ち、sm2z(ti−1) → sgn(sm2z(ti))・|sm2z(ti−1)|とする（処理112）。サイトm3以降も同様な処理を繰り返してN個のsjz(ti)を求め（処理113）、時刻tiの処理が完了する。
【００４９】
時刻tiの処理が完了すれば時刻ti+1の処理に進み、時刻t = tm ≦τまで同様な処理を繰り返して最終的な解を得る。最終的な解は実施例１の場合と同様にsjz(tm)>0ならばsjzfd=1、sjz(tm)<0ならばsjzfd=−1である（処理201及び202）。
【００５０】
スピン反転は量子力学的にはトンネル現象が起こっているといった捉え方もできる。その捉え方に従えば、本実施例はMultiple tunnelingを考慮したものとも言える。
【００５１】
得られた解のエネルギー値は式(7)
【００５２】
【数７】
2019079225.tif 000009
【００５３】
(7)
により与えられる。局所場応答法は量子アニールと類似の動作をするものであり、t = 0で基底状態に準備したスピン系を時間発展させて理想的にはt =τで問題を設定した系の基底状態に導く。本実施例の方法は解の収束性が高く、t = tmの状態は高い確率で演算中の最低エネルギー状態（高い確率で基底状態）である。しかし、実施例１の方法では解の収束性が悪く、最低エネルギー状態がt < tmで現れることもあり、実施例１を用いた場合、各時刻で式(7)を用いてエネルギーを算出し、その中で最低エネルギーであった状態を最終解に選ぶ必要がある。そのための計算量は式(7)の第1項がO(N2)、第2項がO(N)である（Oはランダウの記号）。両者をまとめればO(N2)である。一方、本実施例の方法を用いれば解の収束性が高いのでエネルギーの計算をする必要がなくO(N2)の計算量を節約できる。
【００５４】
その一方で、本実施例の方法は処理111においてスピンの並べ替えを行う。この処理量は以下のように見積もれる。スピンごとに他のスピンと比較する処理を単純に繰り返せばO(N2)である。これが処理量の上限になる。下限は以下のように見積もれる。すでに降べきの順に並んでいるスピン列を考え，隣同士を比較して並べ替えるとする。並べ替えがなければ隣同士の比較だけなので計算量はO(N)である。一般に激しい並べ替えはめったにないので実際の計算量はO(N)に近くなる。よって本実施例のオーバーヘッドはO(N)程度になり実施例１のオーバーヘッドO(N2)に比べて小さい。
【００５５】
局所場応答法全体の計算量は処理102aにおける有効磁場の計算で決まり、Nサイトすべてに対してO(N)の計算をするのでO(N2)である。よってNが十分に大きい系では本実施例のオーバーヘッドO(N)程度は無視できることになる。
【００５６】
以上述べたように、図２の例では処理102a, 103が全てのサイトで独立かつ並列的に処理されていた。そのために処理速度は速いが量子縺れの考慮が不十分であった。一方、本実施例の図４の例では処理112, 111を加えることにより、同時刻の他のスピンの影響を処理102a, 103において反映させ、量子縺れの効果を現象論的に取り入れた。その結果、解精度が向上し、解の収束性が向上し、演算時間が節約できるようになった。
【実施例３】
【００５７】
量子力学的には有効磁場は式(4)に基づき定まる。σ^kzの固有値は±1である。しかし局所場応答法ではスピン変数skzが期待値<σ^kz>の値を取るように動作させるので|skz|≦1である。従って、一般にgjに比べてΣk（≠j）Jkjskzの項を過小評価することになる。
【００５８】
Σk（≠j）Jkjskzの項を過小評価したまま演算させると解精度が下がる。そこでskzの値を参考にgjの値を規格化することにする。gjに因子ci=(Σkskz(ti−1)2/N)1/2を掛けてgjnorm(ti) = cigjとしてgjnorm(ti)を局所項にすればgjnorm(ti)とΣk（≠j）Jkjskzの項の寄与が概ね同等になり、解精度が向上する。尚、離散的に扱う時間軸の分割数をm（tm≦τ）としてc1=1/m程度とする。これはci=(Σkskz(ti−1)2/N)1/2に基づけばskz(t0)=0によりc1=0となってしまうことに対処するためである。
【００５９】
図５に以上の取り扱いを含めたフローチャートを示す。図４との違いは手順102aが手順102bに変わることである。手順102bでは因子ci=(Σkskz(ti−1)2/N)1/2の取り扱いが加わっている。
【００６０】
図６には変形例を示している。パラメタcaを導入し、因子ciをca・ciとすれば因子ciの大きさを調整できる。その場合を図６に示す。後に説明する図９の処理10cも因子ciをca・ciとしている。caの値は経験的に決めるものであり1〜50程度である。
【実施例４】
【００６１】
量子力学的スピン系では常にスピン同士で影響し合っている。即ち、あるサイトjのスピンσ^jzは別のサイトkのスピンσ^kzに影響し、逆にσ^kzがσ^jzに影響する。従って、スピンσ^jzはサイトkのスピンσ^kzを経由して自身に影響する。量子力学においてあるスピンの状態が相互作用相手のスピンの状態だけでなく自身のスピンの状態に依存するのはそのためである。相互作用を通した自身への影響の大きさはΣk（≠j）Jkj2に比例する。実施例１では平均化された量子効果について述べた。Σk（≠j）Jkj2は二乗の項なので平均しても残る。平均化した応答関数がrb0(t)≠1となったのはこのためである。図３の円形の記号及び実線はその様子を示す。Σk（≠j）Jkj2に関する詳細な理論は非特許文献２に記されている。
【００６２】
rb0(t)≠1の原因となったΣk（≠j）Jkj2には問題ごとの情報Jkjが含まれている。この情報を利用できれば演算がより正確になる。そこで平均的応答関数rb0(t)をrb0mod(t)に置き換える。rb0mod(t)は1−rb0mod(t)=(1−rb0(t))Σk（≠j）Jkj2/Σk（≠j）ave(Jkj2)で定義される。ここでave(Jkj2)はrb0(t)を決定する際に使用した問題のJkj2を平均したものである。この改良により応答関数が実際の問題を反映したものになり解精度が向上する。
【００６３】
図７に以上の取り扱いを含めたフローチャートを示す。図５との違いは手順103が手順103cに変わることである。
【実施例５】
【００６４】
量子力学の特徴的な性質として量子縺れと線形重ね合わせが挙げられる。前者の量子縺れの効果は実施例４までに現象論的に取り入れられた。後者の線形重ね合わせの効果は非特許文献２において現象論的に取り入れられている。本実施例では両者を同時に取り入れる。その場合を示したのが図８である。
【００６５】
図８に示す実施例において、図７との違いは手順102bが手順102dに変わることである。
【００６６】
線形重ね合わせが顕著になるのはスピンの符号が変わる時間帯である。その近傍では量子力学的にはバンドが反交差となり、線形重ね合わせ状態となってsjz(t)≒0であり、またそれに連動してBeff,jz(t)≒0になる。この効果を現象論的に取り入れるためには時刻tiにおける有効磁場Beff,jz(ti)を、t=tiとそのひとつ前のt=ti−1における有効磁場の線形結合にすればよい。具体的には時刻ti−1におけるスピンの値sjz(ti−1)により式(8)を使って
【００６７】
【数８】
2019079225.tif 000010
【００６８】
(8)
Bjz0(ti)を求める。次にひとつ前の時刻ti−1における因子も考慮して式(9)により
【００６９】
【数９】
2019079225.tif 000011
【００７０】
(9)
Bjz(ti)を求める。ここでuは0≦u≦1で解精度が高くなるように適当に定める。典型的な値はu≒0.1である。横磁場及びそのスケジュールも含めて有効磁場を記述すれば式(10)になる。
【００７１】
【数１０】
2019079225.tif 000012
【００７２】
(10)
有効磁場が求まれば応答関数rb0mod(t)を利用して
【００７３】
【数１１】
2019079225.tif 000013
【００７４】
(11)
式(11)からsjz(t)を求める。
【実施例６】
【００７５】
実施例５まで、量子縺れの効果を中心に各種量子効果を加えて解精度を向上させる実施例を述べた。しかし、以上の方法を用いても必ず正解に辿り着くとは限らない。そこで本実施例では補助的な手段を述べる。
【００７６】
局所場応答法は決定論的な方法なので同じスピンの初期値で同じ処理をすれば結果はいつも同じである。初期値や処理過程を変更すれば結果が異なる場合もある。そこで初期値や処理過程を変えて複数回磁場挿引の処理を実施し、その中で最低エネルギーを与えた場合のスピン状態を最終解に選んで解精度をさらに向上させる。その場合のフローチャートが図９である。
【００７７】
図９において、処理10a - 10dは図２や図４〜８で示した処理である。処理10a - 10dで得られたスピン値sjzfdからそれぞれエネルギーHpq(tm) (q = 1, 2, ・・・)を処理303で求め、その最低値Ebestを処理304で求める。最低エネルギーを与えたスピン値が最終的な解になる。
【００７８】
処理10aは図２や図４〜８で示した処理をそのまま実行する。処理10bでは処理10aで得られた結果sjzfdの符号反転したものを初期値とする。正確にはパラメタdivを用いて初期値をsjz(t0) = sjzfd/divにする。divで割るのは初期値を十分に小さくするためで、divはm程度の大きさにする。符号反転するのは符号反転した状態がスピンの配位空間上、反転前の状態から最も遠い点であり、局所解に陥った場合に抜け出す機会を与えることができるからである。処理10cでは10aと同様にスピンの初期値を0にするが、局所項gjに掛かる係数ci（実施例３参照）に因子caを掛けてca・ciとすることにより相互作用項と局所項のバランスを変更する。caの値は経験的に決めるものであり1〜50程度である。処理10dはスピンの初期値を乱数で決める。乱数を使うことにより多くの可能性を探れると共に処理を単純化できる。処理10dは乱数の初期値を変えて繰り返し実行する。実行回数が多いほど解精度は向上する。
【実施例７】
【００７９】
本実施例はアルゴリズムとして示されており、通常のコンピュータ上でソフトウェアとして動作させることも専用のハードウェア上で動作させることもできる。本実施例の局所場応答法は比較的単純な演算を繰り返すのが特徴である。従って、繰り返し演算の部分を専用ハードウェアで構築し、その他の部分を汎用装置で実現するのが効果的である。
【００８０】
図１０に本実施例の計算機構成の一例を示す。図１０は通常の計算機の構成と類似であるが局所場応答演算装置600を含む。局所場応答演算装置600は実施例１〜６で述べた演算を専門に行う部分であり、その他の一般的演算は一般演算装置502で行う。
【００８１】
以上の構成は、単体のコンピュータで構成してもよいし、あるいは、主記憶装置501、一般演算装置502、制御装置503、補助記憶装置504、入力装置505、出力装置506等の任意の部分がネットワークで接続された他のコンピュータで構成してもよい。
【００８２】
一般的な演算は通常の計算機と同様な手順で動作させる。記憶部である主記憶装置501と演算部である一般演算装置502間でデータをやり取りし、その繰り返しで演算を進める。その際の司令塔が制御部としての制御装置503である。一般演算装置502で実行されるプログラムは記憶部である主記憶装置501に記憶させる。主記憶装置501で記憶容量が足りない場合は、同じく記憶部である補助記憶装置504を利用する。データやプログラム等の入力には入力装置505を使用し、結果の出力には出力装置506を利用する。入力装置505はキーボードのような手入力装置の他、ネットワーク接続のためのインターフェースも含む。また、このインターフェースは出力装置も兼ねる。
【００８３】
本実施例の局所場応答演算では、実施例１〜６で述べたようにN個のスピン変数sjz(t)とN個の有効磁場変数Beff,jz(t)を交互に繰り返し求める。この繰り返し演算を専門的に実施するのが局所場応答演算装置600である。
【００８４】
実施例６では類似の処理10a - 10dを実施し、処理ごとに解sjzfdを得た。得られた解は局所場応答演算装置600から主記憶装置501にデータ転送し、一般演算装置502を利用してエネルギーHpq(tm)やEbestの計算を行う。即ち、繰り返し演算に属さない個別の演算は一般演算装置502を利用して局所場応答演算装置600の専用性を高める。
【実施例８】
【００８５】
組合せ最適化問題の中で困難度の高い問題はNP困難に属する。またPに分類される問題やNPに分類される問題はすべてNP困難問題に帰着できる。よって、NP困難な組合せ最適化問題を解ければほぼすべての組合せ最適化問題を解けることになる。式(1)の基底状態探索問題はNP困難問題にも対応可能である。本実施例ではその対応の様子を代表的なNP困難問題である最大カット問題を例にして示す。
【００８６】
最大カット問題とはグラフ理論の問題である。グラフ理論ではグラフGを頂点集合Vと辺集合Eから構成しG = (V, E)と書く。辺eは２つの頂点を利用してe = {i, j}と書く。辺eに向きを含めて定義するグラフを有向グラフ、向きの定義を含めないグラフを無向グラフと言う。辺eには重みも定義されておりそれをwij、 wjiと書く。無向グラフならばwij = wjiである。MAX-CUT問題とは、重み付き無向グラフG = (V,  E)の頂点を２つのグループに分ける問題において、カットされる辺の重みの総和を最大化する分割法を求める問題である。分割後の２つの無向グラフをG1 = (V1, E1)及びG2 = (V2,  E2)とすればMAX-CUT問題は
【００８７】
【数１２】
2019079225.tif 000014
【００８８】
を最大化する問題である。頂点i∈V1に対してsi = 1、頂点j∈V2に対してsj = −1とすれば
【００８９】
【数１３】
2019079225.tif 000015
【００９０】
と書ける。最右辺第１項はグラフGが定まれば定数なのでMAX-CUT問題はΣi>jwijsisjを最小化する問題となる。イジング・スピングラスのハミルトニアンは
【００９１】
【数１４】
2019079225.tif 000016
【００９２】
(1)
で与えられる。よってMAX-CUT問題はJij = −wij、 gj = 0とした式(1)の基底状態探索問題と等価になる。
【実施例９】
【００９３】
図４に示した実施例２では、スピンの数であるNサイトの変数全てについて量子縺れを考慮した計算を行っている。しかし、|sjz|の値が大きくなるとスピンが反転する確率が低くなるので、所定以降のスピンについては変化がないものとして量子縺れの効果を無視してもよい。すなわち、図４の処理113をl=x<Nの条件に変更して、Yesの場合にはx+1番目以降のスピンについては処理112を行わずに、処理102a, 103を全てのサイトで独立かつ並列的に処理して時刻tiにおける演算をするように変更しても良い。このようにすると、ある程度の精度を確保しつつ、処理時間を短縮することが可能である。
【００９４】
本実施例中、ソフトウェアで構成した機能と同等の機能は、FPGA（Field Programmable Gate Array）、ASIC（Application Specific Integrated Circuit）などのハードウェアでも実現できる。
【００９５】
本発明は上記した実施形態に限定されるものではなく、様々な変形例が含まれる。例えば、ある実施例の構成の一部を他の実施例の構成に置き換えることが可能であり、また、ある実施例の構成に他の実施例の構成を加えることが可能である。また、各実施例の構成の一部について他の実施例の構成の一部を追加・置換をすると共に、構成の一部を削除ことが可能である。
【符号の説明】
【００９６】
10 - 304 フローチャートにおける各手順を表す
501 主記憶装置
502 一般演算装置
503 制御装置
504 補助記憶装置
505 入力装置
506 出力装置
600 局所場応答演算装置
【図１】
2019079225.tif 000017
【図２】
2019079225.tif 000018
【図３】
2019079225.tif 000019
【図１０】
2019079225.tif 000020
【図４】
2019079225.tif 000021
【図５】
2019079225.tif 000022
【図６】
2019079225.tif 000023
【図７】
2019079225.tif 000024
【図８】
2019079225.tif 000025
【図９】
2019079225.tif 000026