JP2020080006A - シミュレーション装置、コンピュータプログラム及びシミュレーション方法 - Google Patents

Abstract

【課題】最適化問題の近似的解を高速で求めることができるシミュレーション装置、コンピュータプログラム及びシミュレーション方法を提供する。【解決手段】シミュレーション装置は、パウリ行列のｚ成分をスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて目的ハミルトニアンに対する第１確率分布関数を演算する第１演算部と、パウリ行列のｘ成分をトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて初期ハミルトニアンに対する第２確率分布関数を演算する第２演算部と、第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより目的ハミルトニアンのスピン配位の期待値を算出する期待値算出部とを備える。【選択図】図１

Description

本発明は、シミュレーション装置、コンピュータプログラム及びシミュレーション方法に関する。

最適化問題に対して、量子アニーリングと呼ばれる汎用的解法が存在する。量子アニーリングでは、解きたい最適化問題に対して、重ね合わせの状態を生み出す量子揺らぎを印加する。さらに、量子揺らぎの導入が可能な量子コンピュータが開発されているが、扱えるビット数が、例えば２０４８量子ビットであり、スピン数は２０４８個までに限られ、また２体相互作用しか考慮されていないので、実用上の制限が大きい。このため、大規模な実問題に対しては、通常のデジタルコンピュータ上でのシミュレーションが鍵となる。

特許文献１には、従来、デジタルコンピュータ上ではシミュレーションすることができなかった、反強磁性ＸＸ相互作用、いわゆる横磁場以外の量子力学的効果が設定されたシミュレーションをデジタルコンピュータ上で実施できる適応的量子モンテカルロ法が開示されている。

特開２０１８−０６７２００号公報

しかし、特許文献１の方法は、精度の良い解は得られるものの、大規模な最適化問題では、計算時間が長くなる。このため、精度は多少犠牲にしても短時間で計算結果が得られることが望まれる。

本発明は斯かる事情に鑑みてなされたものであり、最適化問題の近似的解を高速で求めることができるシミュレーション装置、コンピュータプログラム及びシミュレーション方法を提供することを目的とする。

本発明の実施の形態に係るシミュレーション装置は、イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートするシミュレーション装置であって、前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算する第１演算部と、前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する第２演算部と、前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する期待値算出部とを備える。

本発明の実施の形態に係るコンピュータプログラムは、コンピュータに、イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算する処理と、前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する処理と、前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する処理とを実行させる。

本発明の実施の形態に係るシミュレーション方法は、イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートするシミュレーション方法であって、前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算し、前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算し、前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する。

本発明によれば、大規模な最適化問題の近似的解を高速に算出することができる。

本実施の形態のシミュレーション装置の構成の一例を示す説明図である。 イジング模型のスピン間の相互作用の一例を示す模式図である。 鈴木トロッター分解の一例を示す模式図である。 相互作用ノードでのメッセージの処理の一例を示す模式図である。 変数ノードでのメッセージの処理の一例を示す模式図である。 確率伝播アルゴリズムに用いる拡張されたグラフ構造の一例を示す模式図である。 拡張されたグラフ構造の相互作用ノードでのメッセージの処理の一例を示す模式図である。 拡張されたグラフ構造の変数ノードでのメッセージの処理の一例を示す模式図である。 本実施の形態のシミュレーション装置による量子アニーリング処理手順の一例を示すフローチャートである。 本実施の形態のシミュレーション装置の構成の他の例を示す説明図である。

以下、本発明をその実施の形態を示す図面に基づいて説明する。図１は本実施の形態のシミュレーション装置１００の構成の一例を示す説明図である。シミュレーション装置１００は、装置全体を制御する制御部１０、入力部１１、目的ハミルトニアン演算部１２、初期ハミルトニアン演算部１３、密度行列算出部１４、確率分布関数演算部１５、出力部１８、確率伝播処理部１９、期待値算出部２０、磁化変数算出部２１、及び記憶部２２を備える。確率分布関数演算部１５は、第１演算部１６及び第２演算部１７を備える。

入力部１１は、シミュレーションを実行するための入力データ及びパラメータを取得する。入力データは、例えば、最適化問題（「組み合せ最適化問題」ともいう）をイジング模型に表現（翻訳）したデータである。イジング模型は、磁性体の振る舞いを記述する数理的な模型であり、±１という２つの値をとるスピン変数σ_i (＝±１)がサイトと称される格子点に配置され、この変数σ_i は磁気モーメントの向きを示す微視的変数として利用される。パラメータは、例えば、温度、横磁場による量子力学的効果、横磁場以外の量子力学的効果、トロッター数などを含む。

出力部１８は、シミュレーションの結果である出力データを出力する。出力データは、最適化問題が表現されたイジング模型の各スピン変数の期待値（スピン配位の期待値）である。

以下、目的ハミルトニアンの算出方法について説明する。

式（１）は、最適化問題をイジング模型で表現したものであり、コスト関数とも称される。

式（１）において、σ_i は二値変数（例えば、±１）であり、添え字のｉはイジング模型の自由度であるスピンが配置された場所を表す。μはスピン間の相互作用を表す添え字であり、Ｊ_μは相互作用の強さを表す。また、ｉ∈∂μは、「μの周りのｉ」という意味である。式（１）の意味を簡単な例で説明する。

図２はイジング模型のスピン間の相互作用の一例を示す模式図である。図２に示すように、スピンをσ₁ からσ₈ で表す（スピン数８）。σ₁ からσ₄ の相互作用をμとし、σ₄ とσ₅ の相互作用をμ′とし、σ₄ 及びσ₆ からσ₈ の相互作用をμ′′とすると、イジング模型のハミルトニアン（エネルギー）は、−Ｊ_μσ₁ σ₂ σ₃ σ₄−Ｊ_μ′σ₄ σ₅−Ｊ_μ′′σ₄ σ₆σ₇σ₈で表すことができ、このコスト関数を最小（または最大）にすることにより最適化問題を解くことができる。図２からも分かるように、本実施の形態では、２体相互作用（例えば、Ｊ_μσ₁ σ₂のように、２つのスピンの積で表すことができる）に限定されず、３体以上の相互作用を扱うことができる。また、スピン数Ｎもデジタルコンピュータのメモリ容量を増加すれば、数千以上のスピン数を扱うことができ、大規模な最適化問題を解くことができる。

量子アニーリングでは、最適化問題に対して、重ね合わせの状態を生み出す量子揺らぎを印加する。重ね合わせの状態とは、２つの異なる状態１、２（例えば、スピンσ＝＋１と、σ＝−１）が同時に存在するという量子力学的な状態であり、例えば、繰り返し観測を行った場合に、ある観測時には状態1が観測され、別の観測時には状態２が観測され、状態１、２がそれぞれある確率で存在するという意味である。量子揺らぎを印加するために、以下のような変換処理を行う。すなわち、イジング模型のスピンの二値変数σ_i （＝±１）を、式（２）で示すパウリ行列のｚ成分に対応付ける。

目的ハミルトニアン演算部１２は、入力された最適化問題に基づいて、目的ハミルトニアンＨ₀ ハットを算出する。具体的には、目的ハミルトニアン演算部１２は、式（１）において、式（２）を用いて、式（３）で表される目的ハミルトニアンＨ₀ ハットを算出する。なお、ハット（∧）は行列であることを示す。

初期ハミルトニアン演算部１３は、式（４）の右辺第２項で示すような初期ハミルトニアンを算出する。初期ハミルトニアンは、式（４）の右辺第２項のように、磁場関数ｇで表すことができ、磁場関数ｇの変数には、複数のスピンσ_i （＝±１）に対応させて、式（５）に示すパウリ行列のｘ成分を含む。パウリ行列のｘ成分は、非対角成分を持っているので、ｚ方向に向いたスピンを反転させることができ、量子揺らぎを生み出すことができる。式（４）において、Ｎはスピン数である。

式（６）は、初期ハミルトニアン演算部１３が算出する、初期ハミルトニアンの具体例を示す。式（６）の右辺第１項は、パウリ行列のｘ成分の和の平均成分の一次項を含み、横磁場による量子力学的効果を与えることができる。係数Γは、量子揺らぎの強さを制御するパラメータである。量子アニーリングにおいて、係数Γは、時間の経過とともに小さくなるように制御する。式（６）の右辺第２項は、パウリ行列のｘ成分の和の平均成分の二次項を含み、いわゆる、反強磁性ＸＸ相互作用を表し、横磁場以外の量子力学的効果を与えることができる。γは所定の係数である。なお、磁場関数は、式（６）に限定されるものではない。例えば、磁場関数は、パウリ行列のｘ成分の和の平均成分の２以上の累乗を含めることもできる。

式（４）において、初期状態（時刻ｔ＝０）では、係数Γを非常に大きな値とし、時間の経過とともに係数Γを小さくし、最終的には０にする。最初は大きな量子揺らぎによって数多くの状態の重ね合わせを実現して状態探査をする。各時刻における瞬間的な基底状態を連続的にたどり、次第にΓが小さくなると、初期ハミルトニアンに比べて目的ハミルトニアンの相対的な重みが大きくなり、最終的には、目的ハミルトニアンの基底状態に到達する。この状態で、最適化問題の解が得られ、スピン配位の期待値を算出することができる。

密度行列算出部１４は、式（４）で示すハミルトニアンＨハットについて、密度行列ρハットを算出する。密度行列ρハットは、式（７）で表すことができる。微視的変数が従う確率分布関数は、ギブスボルツマン分布と呼ばれるが、量子力学では、確率分布の代わりに行列に置き換えられた密度行列が用いられる。

式（７）において、βは、式（８）で表すように、温度Ｔの逆数であり、Ｚは式（９）で表される規格化定数であり、分配関数と呼ばれる。式（９）において、Ｔｒは行列の対角成分の和（トレース）を表す記号である。式（７）のＨハットに式（４）を代入すると式（１０）が得られる。

密度行列算出部１４は、初期ハミルトニアンのパウリ行列のｘ成分の和の平均成分とｘ方向の磁化変数ｍ_x との差を変数とするデルタ関数δと、初期ハミルトニアンを含む指数関数演算子との積を用いて、初期ハミルトニアンに対する密度行列を算出する。具体的には、式（１１）の左辺の初期ハミルトニアンの密度行列は、デルタ関数を用いることにより、式（１１）の右辺のように表すことができる。磁化変数ｍ_x は、イジング模型においてスピンが全体としてどれだけ揃っているかを示す物理量である。

式（１２）はデルタ関数の公式であり、式（１２）を用いることにより、式（１１）の右辺（初期ハミルトニアンの密度行列）は、式（１３）のように表すことができる。すなわち、密度行列算出部１４は、式（１３）で表す、初期ハミルトニアンの密度行列を算出する。式（１２）において、なお、ｘチルダは、物理量ｘに対する演算子であることを示す。上述のように、フーリエ積分を利用して、式（１３）で示すように、初期ハミルトニアンに対する密度行列が、パウリ行列のｘ成分の和の一次項（すなわち、二次項以上の高次項を含まない）だけで表現することができ、横磁場以外の量子力学的効果（反強磁性ＸＸ相互作用など）を横磁場による量子力学的効果のみに置き換えることができる。

式（１４）は、鈴木トロッター分解したときのハミルトニアンＨハット（目的ハミルトニアンＨ₀ ハットと初期ハミルトニアンとの和）に対する密度行列を示す。式（１４）において、τはトロッター数を示す。

図３は鈴木トロッター分解の一例を示す模式図である。図３において、横軸は各サイトに配置されたスピン変数を表し、いわゆる実空間方向を示す。縦軸は鈴木トロッター分解によって導入された方向（トロッター方向）であり、２次元の格子点上に状態変数が配置される。例えば、スピン変数σ_i に対して、トロッター方向に向かって、σ_i1、…、σ_ik、σ_i(k+1)、…、σ_iτが配置される。このように、鈴木トロッター分解によって量子モデルは、次元が一つ増えた状態空間を持つ古典モデルに変換されたと考えることができる。

鈴木トロッター分解により、パウリ行列のｚ成分は、式（１５）のように変換することができ、パウリ行列のｘ成分は、式（１６）のように変換することができる。式（１５）及び式（１６）を用いることにより、式（１４）で表した密度行列は、式（１７）のように表すことができる。すなわち、密度行列算出部１４は、式（１７）で表す密度行列を算出する。テイラー展開、代数的対応関係を用いてイジング模型を仮想的な虚時間方向（トロッター方向）へ拡張して、最適化問題を書き換えることができる。式（１７）において、αは、後述の確率伝播アルゴリズムで繰り返し処理を行う際の更新変数であり、更新変数αは、式（１８）のように更新関数{ tanh（β・ｍ_xチルダ／τ）}に基づいて、更新関数の関数値として求めることができる。

次に、確率伝播アルゴリズムについて説明する。確率伝播法（ＢＰ：Belief Propagation、誤差伝播法ともいう）は、複数の確率変数と密度関数との依存関係をノードで接続したグラフで記述し、そのグラフ構造を利用して高速に確率分布を推論するものであり、そのグラフ上での局所的なメッセージの交換及び処理を行うことにより大域的（全体）の確率分布を推論する。

図４は相互作用ノードでのメッセージの処理の一例を示す模式図である。図４では、便宜上、４つの変数ノードと、４つの変数ノードに繋がる相互作用ノードμとを図示している。相互作用ノードμから変数ノードｉへのメッセージＭμ→ｉは、式（１９）で表すことができる。

式（１９）において、∂μ／ｉは、「ｉを除くμの周り」という意味であり、ｌ∈∂μ／ｉは、図４の例では、変数ノードｌ１、ｌ２、ｌ３を示す。図４に示すように、メッセージＭμ→ｉは、変数ノードｌ（＝ｌ１、ｌ２、ｌ３）から相互作用ノードμへのメッセージＭｌ→μの和に、相互作用ノードμの密度関数ｆｕを積算して求めることができる。

密度関数ｆｕは、式（２０）で表すことができ、式（１）から分かるように、最適化問題のコスト関数が表現されている。

図５は変数ノードでのメッセージの処理の一例を示す模式図である。図５では、便宜上、４つの相互作用と、４つの相互作用ノードに繋がる変数ノードｉとを図示している。変数ノードｉから相互作用ノードμへのメッセージＭｉ→μは、式（２１）で表すことができる。

式（２１）において、∂ｉ／μは、「μを除くｉの周り」という意味であり、ν∈∂ｉ／μは、図５の例では、相互作用ノードν１、ν２、ν３を示す。図５に示すように、メッセージＭｉ→μは、相互作用ノードν（＝ν１、ν２、ν３）から変数ノードｉへのメッセージＭν→ｉの和に、変数ノードｉの密度関数ｆｉを積算して求めることができる。

密度関数ｆｉは、式（２２）で表すことができ、変数ノードｉでは、特に何の処理もせずに、メッセージＭν→ｉの和をメッセージＭｉ→μとして相互作用ノードμへ伝播させる。変数ノードの各変数に適当な初期値を与えて、メッセージの算出及び伝播を繰り返すことにより、メッセージの値が収束し、コスト関数を最小（又は最大）にする変数を求めることができる。

次に、鈴木トロッター分解によって書き換えられた最適化問題に対して確率伝播アルゴリズム（拡張レベルの確率伝播アルゴリズム）を適用する方法について説明する。

図６は確率伝播アルゴリズムに用いる拡張されたグラフ構造の一例を示す模式図である。図６の左側の図は、１つの相互作用ノードμの周りに４つの変数ノードが接続されたグラフ構造を示す。４つの変数ノードのうち、１つの変数ノードをσ_i で表す。変数ノードσ_i には、式（３）で表した、目的ハミルトニアンＨ₀ハットに含まれるσ_iのｚ成分と、式（６）で表した、初期ハミルトニアンに含まれるσ_iのｘ成分との行列成分が存在する。

図６の右側の図は、鈴木トロッター分解により、左側のグラフ構造をトロッター方向（虚時間方向）に拡張したものである。１つの相互作用ノードμの周りに４つの変数ノードが接続されたグラフ構造をトロッター数τだけコピーしたものである。ここで、ｋ＝１、…、ｋ、ｋ＋１、…、ｋ＝τ（トロッター数τ）としている。鈴木トロッター分解により拡張されたグラフ構造では、変数ノードσ_i には、σ_ikのトロッター数τ分の和（σ_iのｚ成分を置き換えたもの）、及びσ_ikσ_i(k+1) のトロッター数τ分の和（σ_iのx成分を置き換えたもの）の変数が存在する。ここで、σ_ikσ_i(k+1)は、ｋ番目の変数ノードのσ_i と（ｋ＋１）番目の変数ノードのσ_i(k+1) との相互作用を表す。

図７は拡張されたグラフ構造の相互作用ノードでのメッセージの処理の一例を示す模式図である。図７では、便宜上、相互作用ノードμと、変数ノードｉとのグラフ構造を、ｋ＝１、２、３、…τだけ拡張したものを図示している。相互作用ノードμから変数ノードｉへのメッセージＭμ→ｉは、式（２３）で表すことができる。

式（２３）において、密度関数ｆｕは、式（２４）で表すことができる。

式（２４）に示すように、第１演算部１６は、パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数σ_ikに置き換え、置き換えたスピン変数σ_ikが含まれる指数関数演算子を用いて目的ハミルトニアンに対する第１確率分布関数としての密度関数ｆｕを演算する。

式（２４）と式（２０）とを対比すると、式（２０）のσ_i が式（２４）ではσ_ikに置き換わっている。また、図７及び式（２３）に示すように、相互作用ノードμから変数ノードｉへのメッセージＭμ→ｉは、相互作用ノードμへのメッセージＭｌ→μの和と密度関数ｆｕとの積をトロッター数τ分だけ加えたものとなっている。

図８は拡張されたグラフ構造の変数ノードでのメッセージの処理の一例を示す模式図である。図８では、図７と同様に、便宜上、相互作用ノードμと、変数ノードｉとのグラフ構造を、ｋ＝１、２、３、…τだけ拡張したものを図示している。変数ノードｉから相互作用ノードμへのメッセージＭｉ→μは、式（２５）で表すことができる。

式（２５）において、密度関数ｆｉは、式（２６）で表すことができる。

式（２６）に示すように、第２演算部１７は、パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数σ_ik、σ_i(k+1)に置き換え、置き換えたスピン変数の積σ_ikσ_i(k+1) をトロッター方向に亘って合計した値を含む指数関数演算子を用いて初期ハミルトニアンに対する第２確率分布関数としての密度関数ｆｉを演算する。図８及び式（２６）に示すように、変数ノードｉでは、トロッター方向に沿って隣り合うσ_ikσ_i(k+1) の相互作用が考慮されている。

第２演算部１７は、密度行列算出部１４が算出した、式（１１）で表す、密度行列に対して、鈴木トロッター分解を行って第２確率分布関数としての密度関数ｆｉを演算することができる。

また、式（２６）には、更新変数αが含まれている。更新変数αは、式（１８）で表す更新関数{ tanh（β・ｍ_xチルダ／τ）}の関数値とすることができ、式（１８）のｍ_xチルダは、後述の式（２８）に示すように、ｘ方向の磁化変数ｍ_xを変数とする磁場関数ｇ（ｍ_x）の導関数ｇ′（ｍ_x）である。すなわち、第２演算部１７は、ｘ方向の磁化変数を変数とする磁場関数の導関数を変数とする更新関数に基づく更新変数αをさらに含む指数関数演算子を用いて初期ハミルトニアンに対する第２確率分布関数としての密度関数ｆｉを演算する。

確率伝播処理部１９は、式（２３）で表す第１メッセージとしての相互作用ノードでのメッセージＭμ→ｉを生成する。相互作用ノードでのメッセージＭμ→ｉは、密度関数ｆｕをトロッター方向に亘って合計した確率分布関数を含む。また、確率伝播処理部１９は、式（２５）で表す第２メッセージとしての変数ノードでのメッセージＭｉ→μを生成する。変数ノードでのメッセージＭｉ→μは、スピン変数の積σ_ikσ_i(k+1) をトロッター方向に亘って合計した値を含む指数関数演算子を有する密度関数ｆｉを含む。

確率伝播処理部１９は、メッセージＭμ→ｉ、及びＭｉ→μの算出及び伝播を繰り返す処理を行う。

期待値算出部２０は、変数ノードの各スピンに適当な初期値を与えて、確率伝播処理部１９によるメッセージの算出及び伝播の処理を繰り返し、メッセージＭμ→ｉ、及びＭｉ→μの値が収束したときの、最適化問題についてのコスト関数を最小（又は最大）にするスピン変数（スピン配位の期待値）を算出する。

次に、確率伝播アルゴリズムでの更新変数αの決定方法について説明する。

磁化変数算出部２１は、更新関数を底とし、トロッター方向に沿って隣り合うスピン変数の積σ_ikσ_i(k+1) を冪指数とする冪演算をスピン数に亘って合計するとともにトロッター数に亘って合計してｘ方向の磁化変数ｍ_xを算出する。磁化変数ｍ_xは、式（２７）により算出することができる。式（２７）において、tanh（β・ｍ_xチルダ／τ）が更新関数である。

式（１８）により、更新変数αは、更新関数の対数によって得られる関数値で表すことができる。

式（２７）において、ｍ_xチルダは、式（２８）により求めることができる。ｍ_xチルダは、磁化変数ｍ_xを変数とする磁場関数ｇ（ｍ_x）の導関数で表すことができる。

メッセージを算出し伝播させたときのσ_ik、σ_i(k+1)（＋１、又は−１）を式（２７）に代入して、磁化変数ｍ_xを算出する。算出した磁化変数ｍ_xを式（２８）に代入して、ｍ_xチルダを算出する。算出したｍ_xチルダを式（１８）に代入して、更新変数αを算出する。算出した更新変数αを式（２６）に代入して、再度、メッセージＭｉ→μとメッセージＭμ→ｉとを算出し伝播処理を行う。なお、メッセージＭｉ→μが算出できると、式（２３）により、メッセージＭμ→ｉを算出できる。伝播処理によって得られたσ_ik、σ_i(k+1)（＋１、又は−１）を式（２７）に代入して、以下、同様の処理を、メッセージの値が収束するまで繰り返す。メッセージの値が収束すると、スピン配位の期待値の近似的解を得ることができる。

記憶部２２は、入力データ、シミュレーション中に得られた処理結果、出力データなどを記憶することができる。

次に、本実施の形態のシミュレーション装置１００の動作について説明する。図９は本実施の形態のシミュレーション装置１００による量子アニーリング処理手順の一例を示すフローチャートである。以下では、便宜上、処理の主体を制御部１０として説明する。制御部１０は、最適化問題をイジング模型で表現したデータを取得し（Ｓ１１）、パラメータを設定する（Ｓ１２）。パラメータは、温度、量子揺らぎ作用を決定付ける磁場関数を表現するデータであり、横磁場、反強磁性ＸＸ相互作用、パウリ行列のｘ成分の和の３次以上の累乗の項などを含む。また、パラメータには、トロッター数が含まれる。

制御部１０は、目的ハミルトニアンを算出し（Ｓ１３）、初期ハミルトニアンを算出する（Ｓ１４）。目的ハミルトニアンは、式（３）により算出することができ、初期ハミルトニアンは、例えば、式（６）により算出することができる。なお、初期ハミルトニアンは、式（６）で表現されるものに限定されない。初期ハミルトニアンに、パウリ行列のｘ成分の和の平均成分の三次以上の累乗の項を含めることもできる。

制御部１０は、初期ハミルトニアンのパウリ行列のｘ成分の和の平均成分を磁化変数ｍ_xで置き換え（Ｓ１５）、目的ハミルトニアンに対する第１確率分布関数を算出し（Ｓ１６）、初期ハミルトニアンに対する第２確率分布関数を算出する（Ｓ１７）。第１確率分布関数は、式（２４）で表され、第２確率分布関数は、式（２６）で表される。

制御部１０は、第１メッセージＭμ→ｉを算出する（Ｓ１８）。第１メッセージＭμ→ｉは、式（２３）で表される。制御部１０は、磁化変数ｍ_xを算出する（Ｓ１９）。磁化変数ｍ_xは、式（２７）で表される。制御部１０は、更新変数αを算出する（Ｓ２０）。更新変数αは、式（１８）で表される。

制御部１０は、第２メッセージＭｉ→μを算出する（Ｓ２１）。第２メッセージＭｉ→μは、式（２５）で表される。制御部１０は、第１メッセージ及び第２メッセージを用いて確率伝播処理を行い（Ｓ２２）、メッセージの値が収束したときのスピン配位の期待値を算出し（Ｓ２３）、処理を終了する。

図１０は本実施の形態のシミュレーション装置の構成の他の例を示す説明図である。図１０において、符号３００は、通常のコンピュータである。コンピュータ３００は、制御部３０、入力部４０、出力部５０、外部Ｉ／Ｆ（インタフェース）部６０などを備える。制御部３０は、ＣＰＵ３１、ＲＯＭ３２、ＲＡＭ３３、Ｉ／Ｆ（インタフェース）３４などを備える。

入力部４０は、シミュレーションのための入力データを取得する。出力部５０は、シミュレーション結果である出力データを出力する。Ｉ／Ｆ３４は、制御部３０と、入力部４０、出力部５０及び外部Ｉ／Ｆ部６０それぞれとの間のインタフェース機能を有する。

外部Ｉ／Ｆ部６０は、コンピュータプログラムを記録した記録媒体Ｍ（例えば、ＤＶＤなどのメディア）からコンピュータプログラムを読み取ることが可能である。

なお、図示していないが、記録媒体Ｍに記録されたコンピュータプログラムは、持ち運びが自由なメディアに記録されたものに限定されるものではなく、インターネット又は他の通信回線を通じて伝送されるコンピュータプログラムも含めることができる。また、コンピュータには、複数のプロセッサを搭載した１台のコンピュータ、あるいは、通信ネットワークを介して接続された複数台のコンピュータで構成されるコンピュータシステムも含まれる。

本実施の形態のシミュレーション装置及びシミュレーション方法によれば、鈴木トロッター分解による虚時間方向へのイジング模型の拡張を行って最適化問題を書き換え、書き換えられた最適化問題に対して確率伝播アルゴリズムを適用した拡張レベルの確率伝播アルゴリズムを用いるので、量子アニーリングマシン（量子コンピュータ）を介さずに、通常のデジタルコンピュータ上で、大規模な最適化問題（組み合わせ最適化問題）の近似的解を高速に求めることができる。

本実施の形態のシミュレーション装置及びシミュレーション方法によれば、反強磁性ＸＸ相互作用を含む、パウリ行列のｘ成分の和の平均成分の三次以上の累乗の項を含む量子力学的効果（量子揺らぎ作用）を、横磁場のみに置き換えることができるので、本実施の形態のシミュレーション装置及びシミュレーション方法によって得られた大規模な最適化問題の解の候補（第１段階の解）を、横磁場のみによる量子力学的効果を扱うことができる量子コンピュータへ提供して、２段階の解法を行うこともできる。

本実施の形態のシミュレーション装置は、パーソナルコンピュータ、ワークステーション、サーバ、ＧＰＵ、ＦＰＧＡ等で実装することができる。

本実施の形態のシミュレーション装置は、イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートするシミュレーション装置であって、前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算する第１演算部と、前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する第２演算部と、前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する期待値算出部とを備える。

本実施の形態のコンピュータプログラムは、コンピュータに、イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算する処理と、前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する処理と、前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する処理とを実行させる。

本実施の形態のシミュレーション方法は、イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートするシミュレーション方法であって、前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算し、前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算し、前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する。

第１演算部は、パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数σ_i が含まれる指数関数演算子を用いて目的ハミルトニアンに対する第１確率分布関数ｆｕを演算する。

第２演算部は、パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数σ_ik、σ_i(k+1)に置き換え、置き換えたスピン変数の積σ_ikσ_i(k+1)をトロッター方向に亘って合計した値を含む指数関数演算子を用いて初期ハミルトニアンに対する第２確率分布関数ｆｉを演算する。

期待値算出部は、第１確率分布関数ｆｕをトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び第２確率分布関数ｆｉを含む第２メッセージを用いた確率伝播アルゴリズムにより目的ハミルトニアンのスピン配位の期待値を算出する。

スピンが配置されるサイトを変数ノードｉとし、スピン間の相互作用を相互作用ノードμとし、変数ノードｉ及び相互作用ノードμを繋ぐグラフ構造上での確率伝播アルゴリズムにおいて、相互作用ノードμから変数ノードｉへのメッセージを第１メッセージＭμ→ｉとし、変数ノードｉから相互作用ノードμへのメッセージを第２メッセージＭｉ→μとする。

第１メッセージＭμ→ｉは、相互作用ノードμでの第１確率分布関数ｆｕをトロッター方向に亘って合計した確率分布関数を含む。第２メッセージＭｉ→μは、変数ノードｉでの第２確率分布関数ｆｉを含む。

第２確率分布関数ｆｉにより、パウリ行列のｘ成分を持つ初期ハミルトニアンで表す量子力学的効果（量子揺らぎ作用）を、トロッター間相互作用に置き換えることができ、数値計算を実行することが可能となる。また、鈴木トロッター分解によって、仮想的な虚時間方向（トロッター方向）へのイジング模型の拡張を行い、拡張したグラフ構造上での確率伝播アルゴリズム（拡張レベルの確率伝播アルゴリズム）を用いることにより、最適化問題の近似的解を高速で求めることができる。

本実施の形態のシミュレーション装置は、前記初期ハミルトニアンの前記パウリ行列のｘ成分の和の平均成分とｘ方向の磁化変数との差を変数とするデルタ関数と、前記初期ハミルトニアンを含む指数関数演算子との積を用いて前記初期ハミルトニアンに対する密度行列を算出する密度行列算出部を備え、前記第２演算部は、前記密度行列に対して鈴木トロッター分解を行って前記第２確率分布関数を演算する。

密度行列算出部は、初期ハミルトニアンのパウリ行列のｘ成分の和の平均成分とｘ方向の磁化変数との差を変数とするデルタ関数と、初期ハミルトニアンを含む指数関数演算子との積を用いて初期ハミルトニアンに対する密度行列を算出する。

これにより、初期ハミルトニアンに対する密度行列が、パウリ行列のｘ成分の和の一次項（すなわち、二次項以上の高次項を含まない）だけで表現することができ、横磁場以外の量子力学的効果（反強磁性ＸＸ相互作用など）を横磁場による量子力学的効果のみに置き換えることができる。

第２演算部は、密度行列に対して鈴木トロッター分解を行って第２確率分布関数ｆｉを演算する。第２確率分布関数ｆｉにより、パウリ行列のｘ成分を持つ初期ハミルトニアンで表す横磁場以外の量子力学的効果（反強磁性ＸＸ相互作用など）を横磁場による量子力学的効果のみに置き換えることができ、置き換えられた横磁場による量子力学的効果（量子揺らぎ作用）を、トロッター間相互作用に置き換えることができ、数値計算を実行することが可能となる。

本実施の形態のシミュレーション装置において、前記第２演算部は、ｘ方向の磁化変数を変数とする磁場関数の導関数を変数とする更新関数に基づく更新変数をさらに含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する。

第２演算部は、ｘ方向の磁化変数ｍ_xを変数とする磁場関数ｇ（ｍ_x）の導関数ｇ′（ｍ_x）を変数とする更新関数（例えば、{ tanh（β・ｇ′（ｍ_x）／τ）}に基づく更新変数α（例えば、α＝ｌｎ{ tanh（β・ｇ′（ｍ_x）／τ）}をさらに含む指数関数演算子を用いて初期ハミルトニアンに対する第２確率分布関数を演算する。βは温度Ｔの逆数であり、τはトロッター数である。

トロッター間相互作用σ_ikσ_i(k+1)が変わると磁化変数ｍ_xが変わり、磁化変数ｍ_xが変わると導関数ｇ′（ｍ_x）の値が変わり、導関数ｇ′（ｍ_x）の値が変わると、更新変数αが変わる。更新変数αが変わると、第２確率分布関数ｆｉが変わり、第２確率分布関数ｆｉが変わると、第２メッセージＭｉ→μ及び第１メッセージＭμ→ｉが変わる。第２メッセージＭｉ→μ及び第１メッセージＭμ→ｉが変わると、トロッター間相互作用σ_ikσ_i(k+1)が変わる。このような更新を繰り返し、メッセージの値が収束すると、各スピン変数のスピン配位の期待値を算出することができる。

本実施の形態のシミュレーション装置は、前記更新関数を底とし、トロッター方向に沿って隣り合うスピン変数の積を冪指数とする冪演算をスピン数に亘って合計するとともにトロッター数に亘って合計してｘ方向の磁化変数を算出する磁化変数算出部を備える。

磁化変数算出部は、更新関数（例えば、{ tanh（β・ｍ_xチルダ／τ）}を底とし、トロッター方向に沿って隣り合うスピン変数の積σ_ikσ_i(k+1)を冪指数とする冪演算をスピン数（例えば、Ｎ）に亘って合計するとともにトロッター数τに亘って合計してｘ方向の磁化変数ｍ_xを算出する。ここで、ｍ_xチルダ＝ｇ′（ｍ_x）である。ｇ（ｍ_x）は磁場関数である。これにより、トロッター間相互作用σ_ikσ_i(k+1)が変わると磁化変数ｍ_xが変わることを定式化できる。

本実施の形態のシミュレーション装置において、前記磁場関数は、前記パウリ行列のｘ成分の和の平均成分の２以上の累乗の項を含む。

磁場関数ｇ（ｍ_x）は、パウリ行列のｘ成分の和の平均成分の２以上の累乗の項を含む。２乗の場合には、量子力学的効果は、反強磁性ＸＸ相互作用による効果とすることができる。これにより、横磁場以外の量子力学的効果を考慮したシミュレーションを行うことができ、汎用性の高い最適化問題の近似的解を高速で求めることができる。

上述の実施の形態では、変数ノードと相互作用ノードとで構成されるグラフ構造上での拡張された確率伝播アルゴリズムについて説明したが、グラフ構造はこれに限定されるものではなく、例えば、ベイジアンネットワークにおいても拡張された確率伝播アルゴリズムを適用して同様の結果を得ることができる。

１０、３０ 制御部
１１、４０ 入力部
１２ 目的ハミルトニアン演算部
１３ 初期ハミルトニアン演算部
１４ 密度行列算出部
１５ 確率分布関数演算部
１６ 第１演算部
１７ 第２演算部
１８、５０ 出力部
１９ 確率伝播処理部
２０ 期待値算出部
２１ 磁化変数算出部
２２ 記憶部
３１ ＣＰＵ
３２ ＲＯＭ
３３ ＲＡＭ
３４ Ｉ／Ｆ
６０ 外部Ｉ／Ｆ部

Claims (7)

1. イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートするシミュレーション装置であって、
   前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算する第１演算部と、
   前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する第２演算部と、
   前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する期待値算出部と
   を備えるシミュレーション装置。
2. 前記初期ハミルトニアンの前記パウリ行列のｘ成分の和の平均成分とｘ方向の磁化変数との差を変数とするデルタ関数と、前記初期ハミルトニアンを含む指数関数演算子との積を用いて前記初期ハミルトニアンに対する密度行列を算出する密度行列算出部を備え、
   前記第２演算部は、
   前記密度行列に対して鈴木トロッター分解を行って前記第２確率分布関数を演算する請求項１に記載のシミュレーション装置。
3. 前記第２演算部は、
   ｘ方向の磁化変数を変数とする磁場関数の導関数を変数とする更新関数に基づく更新変数をさらに含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する請求項１又は請求項２に記載のシミュレーション装置。
4. 前記更新関数を底とし、トロッター方向に沿って隣り合うスピン変数の積を冪指数とする冪演算をスピン数に亘って合計するとともにトロッター数に亘って合計してｘ方向の磁化変数を算出する磁化変数算出部を備える請求項３に記載のシミュレーション装置。
5. 前記磁場関数は、前記パウリ行列のｘ成分の和の平均成分の２以上の累乗の項を含む請求項３又は請求項４に記載のシミュレーション装置。
6. コンピュータに、イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートさせるコンピュータプログラムであって、
   コンピュータに、
   前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算する処理と、
   前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算する処理と、
   前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出する処理と
   を実行させるコンピュータプログラム。
7. イジング模型の複数のスピンがパウリ行列のｚ成分で表され、最適化問題を表現する目的ハミルトニアンと、前記複数のスピンに対応するパウリ行列のｘ成分を含み量子揺らぎを表現する初期ハミルトニアンとを時間変化を伴う係数で組み合わせたハミルトニアンに対して、時間経過とともに量子揺らぎを小さくして前記目的ハミルトニアンのスピン配位の期待値をシミュレートするシミュレーション方法であって、
   前記パウリ行列のｚ成分を鈴木トロッター分解によってスピン変数に置き換え、置き換えたスピン変数が含まれる指数関数演算子を用いて前記目的ハミルトニアンに対する第１確率分布関数を演算し、
   前記パウリ行列のｘ成分を鈴木トロッター分解によってトロッター方向に沿って隣り合うスピン変数に置き換え、置き換えたスピン変数の積をトロッター方向に亘って合計した値を含む指数関数演算子を用いて前記初期ハミルトニアンに対する第２確率分布関数を演算し、
   前記第１確率分布関数をトロッター方向に亘って合計した確率分布関数を含む第１メッセージ及び前記第２確率分布関数を含む第２メッセージを用いた確率伝播アルゴリズムにより前記目的ハミルトニアンのスピン配位の期待値を算出するシミュレーション方法。

Images