JP7472062B2 - 計算装置、計算方法およびプログラム - Google Patents

Description

本発明の実施形態は、計算装置、計算方法およびプログラムに関する。

組合せ最適化問題とは、複数の組合せの中から目的に最も適した組合せを選ぶ問題である。組合せ最適化問題は、数学的には、「目的関数」と呼ばれる、複数の離散変数を有する関数を最大化させる問題、または、当該関数を最小化させる問題に帰着される。組合せ最適化問題は、金融、物流、交通、設計、製造、生命科学など各種の分野において普遍的な問題であるが、組合せ数が問題サイズの指数関数のオーダーで増える、いわゆる「組合せ爆発」のため、必ず最適解を算出することができるとは限らない。また、最適解に近い近似解を得ることすら難しい場合が多い。

各分野における問題を解決し、社会のイノベーションおよび科学技術の進歩を促進するために、組合せ最適化問題の解を精度良く算出する技術の開発が要求されている。

特開２０１７－７３１０６号公報 特開２０１９－１４５０１０号公報

A. Lucas, "Ising formulations of many NP problems", Frontiers in Physics 2, 5 (2014), DOI: 10.3389/fphy.2014.00005 M. W. Johnson et al., "Quantum annealing with manufactured spins", Nature 473, 194-198, 11 May. 2011 T. Inagaki et al.,"A coherent Ising machine for 2000-node optimization problems", Science 354, 603 (2016). H. Goto,"Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network", Sci. Rep. 6, 21686 (2016). M. Yamaoka et al., H. IEEE J. Solid-State Circuits 51, 303 (2016). S. Tsukamoto et al., FUJITSU Sci. Tech. J. 53, 8 (2017). Hayato Goto, Kosuke Tatsumura, Alexander R. Dixon, "Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems",Science Advances, Vol. 5, no. 4, eaav2372, 19 Apr. 2019 EGOR S. TIUNOV, ALEXANDER E. ULANOV, AND A. I. LVOVSKY, "Annealing by simulating the coherent Ising machine",OPTICS EXPRESS 10291, Vol. 27, no. 7, Apr. 2019 H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, K. Tatsumura, "High-performance combinatorial optimization based on classical mechanics", Science Advances; 7, eabe7953",Feb. 2021

本発明の実施形態は、組合せ最適化問題の解を精度良く算出する計算装置、計算方法およびプログラムを提供する。

実施形態に係る計算装置は、組合せ最適化問題を解く。前記計算装置は、更新部と、出力部とを備える。前記更新部は、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する。前記出力部は、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する。前記複数の要素は、前記組合せ最適化問題の複数の離散変数に対応する。前記第１変数および前記第２変数のそれぞれは、実数により表される。前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値に変更し、前記第２変数を予め定められた第３値に変更し、前記第１変数が、前記第１値より大きい予め定められた第２値より大きい場合、前記第１変数を前記第２値に変更し、前記第２変数を前記第３値に変更し、前記第２変数に、予め定められた演算により算出された加速値を加算する。前記加速値は、直前時刻における対応する前記第２変数に、所定係数と前記単位時間との積を乗じた値である。

シミュレーテッド分岐アルゴリズムの分岐現象を示す図。 改良アルゴリズムの分岐現象を示す図。 最大カット問題を解いた場合におけるカット数の平均値を示す図。 最大カット問題を解いた場合におけるカット数の最大値を示す図。 最大カット問題を解いた場合における成功確率を示す図。 本実施形態に係る計算装置の機能構成図。 更新部の処理の流れの第１例を示す図。 更新部の処理の流れの第２例を示す図。 更新部のブロック構成図。 情報処理システムの構成図。 管理サーバの構成図。 記憶部に記憶されるデータを示す図。 計算サーバの構成図。 ストレージに記憶されるデータを示す図。

以下、実施形態について説明する。

（組合せ最適化問題）
イジング問題を解くために使われる装置の一例として、イジングマシンが挙げられる。イジングマシンは、イジングモデルの基底状態のエネルギーを計算する。これまで、イジングモデルは、主に強磁性体や相転移現象のモデルとして使われることが多かった。しかし、近年、イジングモデルは、組合せ最適化問題を解くためのモデルとしての利用が増えている。式（１）は、イジングモデルのエネルギーを示す。

ｓ_ｉ、ｓ_ｊはスピンを表す。スピンは、＋１または－１のいずれかの値をとる２値変数である。ｓ_ｉは、ｉ番目のスピンを表す。ｓ_ｊは、ｊ番目のスピンを表す。ｉおよびｊは、１以上、Ｎ以下の任意の整数である。Ｎは、スピンの数を表し、２以上の整数である。ｈ_ｉは、ｉ番目のスピンに作用する局所磁場を表す。Ｊは、２つのスピン間に作用する力を表す結合係数の行列である。Ｊは、対角成分が０である実対称行列である。Ｊ_ｉｊは、Ｊのｉ行ｊ列の要素を表す。つまり、Ｊ_ｉｊは、ｉ番目のスピンと、ｊ番目のスピンとの間に作用する力を表す結合係数である。なお、式（１）のイジングモデルは、スピンの２次式である。イジングモデルは、スピンの３次以上の項を含む拡張されたモデル（多体相互作用を有するイジングモデル）であってもよい。

イジングマシンは、式（１）により表されるエネルギーＥ_{Ｉｓｉｎｇ}を目的関数とし、エネルギーＥ_{Ｉｓｉｎｇ}を可能な限り小さくする解を算出する。エネルギーＥ_{Ｉｓｉｎｇ}が最小値となるイジングモデルの解（ｓ_１、ｓ_２、・・・、ｓ_Ｎ）は、最適解と呼ばれる。ただし、イジングモデルの解は、最適解でなく、エネルギーＥ_{Ｉｓｉｎｇ}が最小値に近い近似解であってもよい。すなわち、イジング問題は、最適解のみならず、近似解を算出する問題であってもよい。

また、０または１のいずれかの値をとる離散変数（ビット）の２次関数を目的関数とする組合せ最適化問題は、ＱＵＢＯ（Ｑｕａｄｒａｔｉｃ Ｕｎｃｏｎｓｔｒａｉｎｅｄ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）問題と呼ばれる。離散変数（ビット）は、（１＋ｓ_ｉ）／２の演算を用いることにより、ｓ_ｉに変換される。つまり、ＱＵＢＯ問題は、式（１）で表されるイジング問題と等価であるといえる。従って、０－１組合せ最適化問題は、イジング問題に変換し、イジングマシンにより解を算出することが可能である。

イジングマシンは、例えば、量子アニーラ、コヒーレントイジングマシンおよび量子分岐マシン等によりハードウェア実装される。量子アニーラは、超伝導回路を使って量子アニーリングを実現する。コヒーレントイジングマシンは、光パラメトリック発振器で形成されたネットワークの発振現象を利用する。量子分岐マシンは、カー効果を有するパラメトリック発振器のネットワークにおける量子力学的な分岐現象を利用する。これらのハードウェア実装されたイジングマシンは、計算時間の大幅な短縮を実現する可能性がある一方、大規模化および安定的な運用が難しいという課題もある。

イジング問題は、広く普及しているデジタルコンピュータを用いて解を算出することも可能である。デジタルコンピュータは、量子アニーラ、コヒーレントイジングマシンおよび量子分岐マシン等と比べ、大規模化および安定運用が可能である。シミュレーテッドアニーリング（ＳＡ）は、デジタルコンピュータでイジング問題を解くためのアルゴリズムの一例である。ただし、シミュレーテッドアニーリングは、それぞれの変数が逐次更新される逐次更新アルゴリズムであるため、並列化による計算処理の高速化は難しい。

非特許文献７には、組合せ最適化問題を解くためのアルゴリズムとして、シミュレーテッド分岐アルゴリズムが提案されている。シミュレーテッド分岐アルゴリズムは、イジングモデルを用いて、デジタルコンピュータによって、規模の大きい組合せ最適化問題を高速に解くことが可能である。シミュレーテッド分岐アルゴリズムは、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、マイクロプロセッサ、ＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、または、これらの組合せの回路等の電子回路によっても、規模の大きい組合せ最適化問題を高速に解くことが可能である。

（シミュレーテッド分岐アルゴリズム）
シミュレーテッド分岐アルゴリズムは、それぞれがＮ個の要素に対応する２つの変数ｘ_ｉおよび変数ｙ_ｉを用いる。変数ｘ_ｉを第１変数、変数ｙ_ｉを第２変数と呼ぶ場合もある。シミュレーテッド分岐アルゴリズムにおいて、Ｎ個の要素のそれぞれは、仮想的な粒子を表す。Ｎ個の要素は、イジング問題のＮ個のスピンに対応する。従って、Ｎ個の要素は、組合せ最適化問題のＮ個の離散変数（ビット）に対応する。変数ｘ_ｉおよび変数ｙ_ｉは、いずれも、実数で表される連続変数である。変数ｘ_ｉは、Ｎ個の粒子のうちのｉ番目の粒子の位置を表す。変数ｙ_ｉは、ｉ番目の粒子の運動量を表す。Ｎは、イジングモデルに含まれるスピンの数を表し、２以上の整数である。ｉは、１以上、Ｎ以下の任意の整数を表し、Ｎ個の要素のそれぞれを特定するインデックスを表す。

シミュレーテッド分岐アルゴリズムは、それぞれＮ個ある２つの変数ｘ_ｉおよび変数ｙ_ｉについて、下記の式（２）の連立常微分方程式を数値的に解く。

Ｈは、下記の式（３）のハミルトニアンである。

係数Ｄは、予め定められた定数であり、離調（ｄｅｔｕｎｉｎｇ）に相当する。係数ｐ（ｔ）は、ポンピング振幅（ｐｕｍｐｉｎｇ ａｍｐｌｉｔｕｄｅ）に相当し、シミュレーテッド分岐アルゴリズムの計算時に更新回数に応じて値が単調増加する。ｔは、時刻を表す変数である。係数ｐ（ｔ）の初期値は０に設定されていてもよい。係数Ｋは、正のカー係数（Ｋｅｒｒ ｃｏｅｆｆｉｃｉｅｎｔ）に相当する。

ｆ_ｉは、外力を表し、下記の式（４）で表される。

式（４）のｚ_ｉは、式（３）の中の小カッコの内の数式を変数ｘ_ｉで偏微分した式である。式（３）の中の小カッコの内の数式は、イジングモデルのエネルギーＥ_{Ｉｓｉｎｇ}に対応する。

ｃは、係数である。ｃは、例えば、計算を実行する前に予め定められる定数であってもよい。例えば、ｃは、Ｊ行列の最大固有値の逆数に近い値に設定してもよい。また、ｃは、例えば、０．５Ｄ√（Ｎ／２ｎ）であってもよい。ここで、ｎは、組合せ最適化問題に係るグラフのエッジ数である。また、α（ｔ）は、ｐ（ｔ）とともに増加する係数である。例えば、α（ｔ）は、√（ｐ（ｔ））であってもよい。

また、シミュレーテッド分岐アルゴリズムは、３次以上の離散変数を含む目的関数を最小化する組合せ最適化問題を解くことも可能である。３次以上の離散変数を含む目的関数を最小化する組合せ最適化問題は、ＨＯＢＯ（Ｈｉｇｈｅｒ Ｏｒｄｅｒ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）問題と呼ばれる。ＨＯＢＯ問題を解く場合、シミュレーテッド分岐アルゴリズムは、高次へ拡張されたイジングモデルを用いる。

高次へ拡張されたイジングモデルにおけるエネルギーは、下記の式（５）により表される。

ここで、Ｊ^（ｎ）は、ｎ階テンソルを表す。Ｊ^（ｎ）は、式（１）の局所磁場を表すｈ_ｉおよび結合係数の行列を表すＪを一般化して表した表記である。例えば、Ｊ^（１）は、ｈ_ｉに相当する。Ｊ^（ｎ）は、複数の下付き添え字に同じ値を含む場合、要素の値は０となる。式（５）は、３次の項までを表しているが、４次以降も同様に定義することができる。式（５）は、多体相互作用を含むイジングモデルのエネルギーに相当する。

シミュレーテッド分岐アルゴリズムは、ＨＯＢＯ問題を解く場合、式（３）のハミルトニアンＨが式（６）に置き換えられる。

また、シミュレーテッド分岐アルゴリズムは、ＨＯＢＯ問題を解く場合、式（４）の外力ｆ_ｉが式（７）に置き換えられる。

式（７）のｚ_ｉは、式（６）の中の小カッコの内の数式を変数ｘ_ｉで偏微分した式である。式（６）の中の小カッコの内の数式は、高次へ拡張されたイジングモデルのエネルギーＥ_{Ｉｓｉｎｇ}に対応する。

そして、シミュレーテッド分岐アルゴリズムは、ｐ（ｔ）の値を初期値（例えば、０）から所定の値まで増加させた後における変数ｘ_ｉの符号に基づき、スピンｓ_ｉの値を算出する。シミュレーテッド分岐アルゴリズムは、例えば、ｘ_ｉ＞０の場合にｓｇｎ（ｘ_ｉ）＝１、ｘ_ｉ＜０の場合にｓｇｎ（ｘ_ｉ）＝－１となる符号関数を用いて、スピンｓ_ｉの値を算出する。

（シミュレーテッド分岐アルゴリズムの演算）
シミュレーテッド分岐アルゴリズムは、ＱＵＢＯ問題を解く場合、シンプレクティック・オイラー法を用いて、式（２）、式（３）および式（４）によって与えられる微分方程式を解く。ＨＯＢＯ問題を解く場合、シミュレーテッド分岐アルゴリズムは、シンプレクティック・オイラー法を用いて、式（２）、式（６）、式（７）によって与えられる微分方程式を解く。

ここで、シンプレクティック・オイラー法を使う場合、式（２）、式（３）および式（４）によって与えられる微分方程式および式（２）、式（６）、式（７）によって与えられる微分方程式は、式（８）に示すような、離散的な漸化式に書き換えられる。

ｔは、時刻を表す。Δｔは、時間ステップ（単位時間、時間刻み幅）を表す。

シミュレーテッド分岐アルゴリズムを実行する場合、デジタルコンピュータまたはＦＰＧＡ等の電子回路は、式（８）のアルゴリズムに基づき、それぞれＮ個ある２つの変数ｘ_ｉおよび変数ｙ_ｉを、初期時刻から時間ステップ毎に順次に、変数ｘ_ｉと変数ｙ_ｉとを交互に更新する。そして、デジタルコンピュータまたはＦＰＧＡ等の電子回路は、終了時刻におけるＮ個の変数ｘ_ｉの値を、符号関数を用いて２値化して、Ｎ個のスピンの値を出力する。

なお、式（８）は、微分方程式との対応関係を示すために、時刻ｔおよび時間ステップΔｔを用いて表されている。ただし、シンプレクティック・オイラー法をデジタルコンピュータまたはＦＰＧＡ等の電子回路で実行する場合、式（８）を演算するためのアルゴリズムは、明示的なパラメータとして時刻ｔおよび時間ステップΔｔを含まなくてよい。例えば、時間ステップΔｔを１とする場合、式（８）を演算するためのアルゴリズムは、時間ステップΔｔを含まなくてよい。例えば、明示的なパラメータとして時刻ｔを含まない場合、式（８）を演算するアルゴリズムは、ｘ_ｉ（ｔ＋Δｔ）をｘ_ｉ（ｔ）の更新後の値として処理を実行する。すなわち、式（８）を演算するアルゴリズムは、式（８）における“ｔ”を更新前の変数を特定するパラメータ、“ｔ＋Δｔ”を更新後の変数を特定するパラメータとして処理を実行する。以降で説明する式（８）を改良したアルゴリズムも同様である。

（改良アルゴリズム）
発明者は、シミュレーテッド分岐アルゴリズムを改良したアルゴリズム（改良アルゴリズム）を発明した。改良アルゴリズムは、シンプレクティック・オイラー法を用いて、式（８）により表される微分方程式を演算する場合において、次の処理をさらに実行する。

すなわち、改良アルゴリズムは、変数ｘ_ｉの更新後、変数ｘ_ｉが予め定められた第１値より小さい場合、変数ｘ_ｉを第１値に変更し、且つ、変数ｙ_ｉを予め定められた第３値に変更する。また、改良アルゴリズムは、変数ｘ_ｉの更新後、変数ｘ_ｉが第１値より大きい第２値より大きい場合、変数ｘ_ｉを第２値に変更し、且つ、変数ｙ_ｉを第３値に変更する。

例えば、第１値は、－１であり、第２値は、＋１であり、第３値は、０である。この場合、改良アルゴリズムは、変数ｘ_ｉの更新後、変数ｘ_ｉが－１より小さい場合、変数ｘ_ｉを－１に変更し、且つ、変数ｙ_ｉを０に変更する。また、改良アルゴリズムは、変数ｘ_ｉの更新後、変数ｘ_ｉが＋１より大きい場合、変数ｘ_ｉを＋１に変更し、且つ、変数ｙ_ｉを０に変更する。

さらに、改良アルゴリズムは、変数ｘ_ｉおよび変数ｙ_ｉを更新した後に、変数ｙ_ｉに、予め定められた演算により算出された加速値を加算する。加速値は、例えば、直前時刻における変数ｙ_ｉ（すなわち、更新前の変数ｙ_ｉ）に、所定係数（ｒ）と単位時間（Δｔ）との積（ｒΔｔ）を乗じた値である。所定係数（ｒ）は、例えば予め定められた定数である。また、所定係数は、更新毎に乱数に基づき定まる値であってもよい。第１値は、－１であり、第２値は、＋１であり、第３値は、０である場合、所定係数は、例えば、０より大きく、１より小さい値である。

また、式（８）の“Ｋｘ^２ _ｉ（ｔ＋Δｔ）”は、計算中に変数ｘ_ｉが発散するのを防止する機能を有する項である。改良アルゴリズムは、更新毎に、変数ｘ_ｉが第１値（例えば－１）以上、第２値（例えば＋１）以下の範囲に制限されるので、計算中に変数ｘ_ｉが発散することはない。従って、改良アルゴリズムは、式（８）により表される微分方程式における“Ｋｘ^２ _ｉ（ｔ＋Δｔ）”の項を削除した微分方程式を演算することができる。

そこで、改良アルゴリズムは、式（８）のアルゴリズムに代わり、下記の式（９）を用いて演算される。

ＱＵＢＯ問題の場合、式（９）のｚ_ｉ（ｔ＋Δｔ）は、下記の式（１０）により表される。

ＨＯＢＯ問題の場合、式（９）のｚ_ｉ（ｔ＋Δｔ）は、下記の式（１１）により表される。

上述の式（９）の改良アルゴリズムは、式（８）と同様、ハミルトン方程式を解くものであり、変数ｙ_ｉが運動量に相当する。そのため、シンプレクティック・オイラー法を使うことにより、時間ステップΔｔとして小さな値を使わなくても、安定的に解を算出することができる。

なお、改良アルゴリズムによる求解が可能な組合せ最適化問題は、イジング問題に限られず、一般的な２値変数の組合せ最適化問題を解くことが可能である。例えば、改良アルゴリズムは、目的関数の変数が、ａ（第１値）と、ａより大きいｂ（第２値）のいずれかをとる２値変数である、組合せ最適化問題に適用することが可能である。また、一定の更新回数の後に目的関数の解を算出する場合、符号関数の代わりに、値域がａまたはｂの２値である関数ｆ（ｘ_ｉ）を使ってもよい。この関数ｆ（ｘ_ｉ）がとる値は、変数ｘ_ｉの値をしきい値ｖ（ａ＜ｖ＜ｂ）と比較した結果に基づいて決まる。例えば、ｘ_ｉ＜ｖであるならば、ｆ（ｘ_ｉ）＝ａとなる。また、ｖ＜ｘ_ｉであるならば、ｆ（ｘ_ｉ）＝ｂとなる。例えば、ｘ_ｉ＝ｖである場合、ｆ（ｘ_ｉ）＝ａまたは、ｆ（ｘ_ｉ）＝ｂとなる。しきい値ｖは、例えば、（ａ＋ｂ）／２であってもよい。

図１は、シミュレーテッド分岐アルゴリズムにより最適化問題を解いた場合における、変数ｘ_ｉの分岐現象を表す図である。シミュレーテッド分岐アルゴリズムにより最適化問題を解いた場合、系のパラメータが変化することに伴い、安定運動状態が１個のみの系から、安定運動状態が２個の系へと遷移する分岐現象が生じる。図１に示すように、分岐現象が進むと、変数ｘ_ｉは、－１または＋１の近傍に集中するが、－１より小さい領域、または、＋１より大きい領域にも広がる。

図２は、改良アルゴリズムにより最適化問題を解いた場合における、変数ｘ_ｉの分岐現象を表す図である。改良アルゴリズムにより最適化問題を解いた場合、シミュレーテッド分岐アルゴリズムと同様に、系のパラメータが変化することに伴い、安定運動状態が１個のみの系から、安定運動状態が２個の系へと遷移する分岐現象が生じる。しかし、改良アルゴリズムは、－１より小さいｘ_ｉを－１に戻し、＋１より大きいｘ_ｉを＋１に戻すので、－１より小さい領域および＋１より大きい領域に広がらない。

図３および図４は、改良アルゴリズムを用いて、最大カット問題Ｋ２０００を１００００回解いた場合における、時間ステップ数に対するカット数の結果を示す図である。図３は、カット数の平均値を表す。図４は、カット数の最大値を表す。また、図５は、改良アルゴリズムを用いて、最大カット問題Ｋ２０００を１００００回解いた場合における、最大カット数が得られた確率である成功確率を表す。

なお、図３、図４および図５は、Ｄを１、Ｋを１、ｃを０．９√（１９９９）、ｒを０．５、Δｔを１．１に設定し、ｐを０から１に線形に増加させて、改良アルゴリズムを実行した結果を表す。図３および図４の破線Ｃｍａｘは、Ｋ２０００で知られている最大カット数である３３３３７を示している。また、図３および図４を参照すると、改良アルゴリズムは、最大カット数に近づいており、精度良く最適化問題を解くことができることがわかった。また、非特許文献９には、シミュレーテッド分岐アルゴリズムを改良した他の方式（ｂＳＢ，ｄＳＢ）が記載されている。図３、図４および図５には、比較例として、ｂＳＢおよびｄＳＢにより同一のパラメータで最大カット問題Ｋ２０００を解いた結果も示している。改良アルゴリズムは、ｂＳＢおよびｄＳＢと比較しても同等以上の結果が得られることがわかった。また、図５に示すように、改良アルゴリズムは、ｂＳＢおよびｄＳＢと比較して高い成功確率が得られることがわかった。

（機能構成）
図６は、本実施形態に係る計算装置１０の機能構成を示す図である。

計算装置１０は、改良アルゴリズムを用いて、最適化問題を解く。計算装置１０は、例えば、コンピュータ等の情報処理装置、ネットワークを介して複数のコンピュータまたはサーバが相互に通信をして構成されるコンピュータシステム、または、複数台のコンピュータが連携して情報処理を実行するＰＣクラスタ等により実現される。また、計算装置１０は、ＣＰＵ、マイクロプロセッサ、ＧＰＵ、ＦＰＧＡまたはＡＳＩＣ、または、これらの組合せの回路等の電子回路によって実現される。

計算装置１０は、機能構成として、入力部１２と、更新部１４と、出力部１６とを備える。

入力部１２は、組合せ最適化問題の目的関数を定義するための情報（例えば、Ｎ、Ｊ、ｈ）、および、改良アルゴリズムを実行するために必要な係数を表す情報（ｘ、Δｔ、Ｔ，ｐ（ｔ）、α（ｔ））を外部装置から受け取り、更新部１４に与える。更新部１４は、改良アルゴリズムを用いて、第１変数（ｘ_ｉ）および第２変数（ｙ_ｉ）が対応付けられた複数の要素のそれぞれについて、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎に順次に、第１変数（ｘ_ｉ）および第２変数（ｙ_ｉ）を交互に更新する。

そして、出力部１６は、終了時刻（ｔ＝Ｔ）における複数の要素のそれぞれの第１変数（ｘ_ｉ）に基づき、組合せ最適化問題の解を出力する。例えば、出力部１６は、終了時刻における複数の要素のそれぞれについて、第１変数（ｘ_ｉ）を予め設定されたしきい値により２値化した離散変数の値を算出する。そして、出力部１６は、算出した複数の離散変数の値を組合せ最適化問題の解として出力する。

（処理フロー）
図７は、更新部１４の処理の流れの第１例を示す図である。更新部１４は、例えば、図７に示す流れで処理を実行する。

まず、Ｓ１０１において、更新部１４は、パラメータを設定する。具体的には、更新部１４は、Ｎ×Ｎ個の結合係数を含む行列であるＪ、および、Ｎ個の局所磁場を表す局所磁場係数を含む配列であるｈを設定する。なお、更新部１４は、ＨＯＢＯ問題を解く場合には、Ｊおよびｈに代えて、Ｎ^ｍ個の作用係数を含むｍ次テンソルであるＪ^（ｍ）を設定する。この場合、（ｍ）は、ＨＯＢＯ問題の目的関数の変数の次数を表す。さらに、更新部１４は、係数であるＤ、係数であるｃ、単位時間を表すΔｔ、終了時刻を表すＴ、関数であるｐ（ｔ）、および、関数であるα（ｔ）を設定する。ｐ（ｔ）およびα（ｔ）は、ｔ＝初期時刻（例えば０）で０、ｔ＝終了時刻（Ｔ）で１となる増加関数である。更新部１４は、Ｊ、ｈを入力部１２からの受け取った情報に応じて設定する。更新部１４は、Ｄ、ｃ、Δｔ、Ｔ、ｐ（ｔ）およびα（ｔ）を、入力部１２から受け取ったパラメータに応じて設定してもよいし、予め決定されており変更できないパラメータを設定してもよい。

続いて、Ｓ１０２において、更新部１４は、変数を初期化する。具体的には、更新部１４は、時刻を表す変数であるｔを初期時刻（例えば、０）に初期化する。さらに、更新部１４は、Ｎ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））のそれぞれおよびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））のそれぞれに、ユーザから受け取った初期値、予め定められた固定値、または、乱数を代入する。

続いて、更新部１４は、Ｓ１０３とＳ１１８との間のループ処理を、ｔがＴより大きくなるまで繰り返す。１回のループ処理において、更新部１４は、対象時刻（ｔ＋Δｔ）におけるＮ個の第１変数（ｘ_１（ｔ＋Δｔ）～ｘ_Ｎ（ｔ＋Δｔ））を、対象時刻（ｔ）におけるＮ個の第２変数（ｙ_１（ｔ＋Δｔ）～ｙ_Ｎ（ｔ＋Δｔ））に基づき算出する。また、１回のループ処理において、更新部１４は、対象時刻（ｔ＋Δｔ）におけるＮ個の第２変数（ｙ_１（ｔ＋Δｔ）～ｙ_Ｎ（ｔ＋Δｔ））を、直前時刻（ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））および直前時刻（ｔ）におけるＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））に基づき算出する。

なお、直前時刻（ｔ）は、対象時刻（ｔ＋Δｔ）より単位時間（Δｔ）前の時刻である。すなわち、更新部１４は、Ｓ１０３とＳ１１８との間のループ処理を繰り返すことにより、Ｎ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎に順次に更新する。

Ｓ１０４において、更新部１４は、直前時刻（ｔ）におけるＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を記憶する。更新部１４は、記憶したＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を後述のＳ１１５で加速値を算出するために用いる。

続いて、更新部１４は、Ｓ１０５とＳ１０９との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。ｉは、１からＮまでの整数であり、Ｎ個の要素のうちの処理対象を表すインデックスである。Ｎ個の要素のそれぞれは、第１変数（ｘ_ｉ（ｔ））および第２変数（ｙ_ｉ（ｔ））が対応付けられる。Ｓ１０５とＳ１０９との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。

Ｓ１０６において、更新部１４は、Ｎ個の要素のそれぞれの直前時刻（ｔ）における第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ））を算出する。ＱＵＢＯ問題の場合、作用係数は、Ｊに含まれる結合係数およびｈに含まれる局所磁場係数である。ＨＯＢＯ問題の場合、作用係数は、Ｊ^（ｎ）に含まれる。

ＱＵＢＯ問題の場合、更新部１４は、式（１２）を算出する。

また、ＨＯＢＯ問題の場合、更新部１４は、式（１３）を算出する。

続いて、Ｓ１０７において、更新部１４は、更新値（ｚ_ｉ（ｔ））に係数（ｃ）を乗算することにより、外力（ｆ_ｉ（ｔ））を算出する。具体的には、更新部１４は、式（１４）を算出する。

続いて、Ｓ１０８において、更新部１４は、対象要素の対象時刻（ｔ＋Δ）における第２変数（ｙ_ｉ（ｔ＋Δ））を、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に、外力（ｆ_ｉ（ｔ））と対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））とに基づく値に単位時間（Δｔ）を乗算した値を、加算することにより、算出する。具体的には、更新部１４は、式（１５）を算出する。

更新部１４は、以上のようなＳ１０５とＳ１０９との間のループ処理をＮ回実行することにより、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、直前時刻（ｔ）における第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））に基づき更新する。更新部１４は、Ｓ１０５とＳ１０９との間のループ処理をＮ回実行した場合、処理をＳ１１０に進める。

続いて、更新部１４は、Ｓ１１０とＳ１１６との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ１１０とＳ１１６との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。

Ｓ１１１において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））に、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））と予め定められた定数（Ｄ）と単位時間（Δｔ）とを乗算した値を加算することにより算出する。具体的には、更新部１４は、式（１６）を算出する。

続いて、Ｓ１１２において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が予め定められた第１値以上、予め定められた第２値以下であるか否かを判断する。本例において、第１値は、－１であり、第２値は、＋１である。例えば、第２値は、連続量である第１変数（ｘ_ｉ（ｔ））の単位量である。例えば、第１値は、第１値の正負を反転した値である。

更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値以上、第２値以下である場合（Ｓ１１２のＹｅｓ）、処理をＳ１１５に進める（Ｓ１１２のＹｅｓ）。対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値より小さいまたは第２値より大きい場合（Ｓ１１２のＮｏ）、処理をＳ１１３に進める。

Ｓ１１３において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））の制約処理をする。具体的には、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を第１値に変更する。また、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第２値より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を第２値に変更する。本例においては、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が－１より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を－１に変更する。また、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が＋１より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を＋１に変更する。

Ｓ１１３に続いて、Ｓ１１４において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））の制約処理をする。具体的には、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、予め定められた第３値に変更する。第３値は、例えば、第１値より大きく、第２値より小さい値である。本例において、第３値は、０である。更新部１４は、Ｓ１１４を終えると、処理をＳ１１５に進める。

Ｓ１１５において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））に、予め定められた演算により算出された加速値を加算する。例えば、加速値は、直前時刻における第２変数（ｙ_ｉ（ｔ））に、所定係数（ｒ）と単位時間（Δｔ）との積（ｒΔｔ）を乗じた値である。所定係数（ｒ）は、予め定められた定数である。例えば、所定係数（ｒ）は、０より大きく１より小さい値である。この場合、更新部１４は、式（１７）に示す演算を実行する。

なお、所定係数（ｒ）は、時刻（ｔ）に依存しない値であれば、他の値であってもよい。例えば、所定係数（ｒ）は、０より大きく１より小さい範囲内における、乱数に基づく値であってもよい。更新部１４は、Ｓ１１５を終えると、処理をＳ１１６に進める。

更新部１４は、以上のようなＳ１１０とＳ１１６との間のループ処理をＮ回実行することにより、次の処理を実行する。すなわち、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））に基づき更新する。続いて、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値より小さい場合、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を第１値に変更する。また、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第２値より大きい場合、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を第２値に変更する。また、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が、第１値より小さいまたは第２値より大きい場合、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を第３値に変更する。さらに、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））に、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に所定係数（ｒ）と単位時間（Δｔ）との積（ｒΔｔ）を乗じた加速値を加算する。

更新部１４は、Ｓ１１０とＳ１１６との間のループ処理をＮ回実行した場合、処理をＳ１１７に進める。

Ｓ１１７において、更新部１４は、直前時刻（ｔ）に単位時間（Δｔ）を加算して、対象時刻（ｔ＋Δｔ）を更新する。Ｓ１１８において、更新部１４は、Ｓ１０４からＳ１１７までの処理を、ｔが終了時刻（Ｔ）を超えるまで繰り返す。そして、更新部１４は、ｔが終了時刻（Ｔ）より大きくなった場合、本フローを終了する。

そして、更新部１４は、Ｎ個の要素のそれぞれについて、終了時刻（ｔ＝Ｔ）における第１変数（ｘ_ｉ（Ｔ））の符号に応じて、対応するスピンの値を算出する。例えば、更新部１４は、終了時刻（ｔ＝Ｔ）における第１変数（ｘ_ｉ（Ｔ））の符号が負である場合、対応するスピンを－１とし、正である場合、対応するスピンを＋１とする。そして、更新部１４は、算出した複数のスピンの値、または、算出した複数のスピンの値を離散変数に変換した値を組合せ最適化問題の解として出力する。

図７に示すフローチャートに従った処理を実行することにより、更新部１４は、Ｎ個の要素のそれぞれについて、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間毎に順次に、第１変数（ｘ_ｉ（ｔ＋Δｔ））および第２変数（ｙ_ｉ（ｔ＋Δｔ））を交互に更新する。また、図７に示すフローチャートに従った処理を実行することにより、更新部１４は、単位時間毎の更新処理において、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を算出した後に、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を算出する。さらに、図７に示すフローチャートに従った処理を実行することにより、更新部１４は、単位時間毎の更新処理において、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を更新する前において、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を算出し、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を更新した後において、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））に加速値を加算する。

以上により、更新部１４は、シンプレクティック・オイラー法を用いて改良アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

図８は、更新部１４の処理の流れの第２例を示す図である。更新部１４は、改良アルゴリズムを用いて最適化問題を解く場合、図７に示す流れに代えて、図８に示す流れで処理を実行してもよい。なお、図８のフローチャートは、図７で示したフローチャートの処理とほぼ同一の処理を実行するステップには、同一のステップ番号を付けている。

まず、Ｓ１０１およびＳ１０２において、更新部１４は、図７に示す第１例と同様の処理を実行する。続いて、更新部１４は、Ｓ１０３とＳ１１８との間のループ処理を、図７に示す第１例と同様に実行する。

Ｓ１０４において、更新部１４は、図７に示す第１例と同様の処理を実行する。

続いて、更新部１４は、Ｓ１１０とＳ１１６との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ１１０とＳ１１６との間のループ処理において、更新部１４は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。

Ｓ１１１において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））に、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））と予め定められた定数（Ｄ）と単位時間（Δｔ）とを乗算した値を加算することにより算出する。具体的には、更新部１４は、式（１８）を算出する。

続いて、Ｓ１１２において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が予め定められた第１値以上、予め定められた第２値以下であるか否かを判断する。更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値以上、第２値以下である場合（Ｓ１１２のＹｅｓ）、処理をＳ１１５に進める（Ｓ１１２のＹｅｓ）。対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値より小さいまたは第２値より大きい場合（Ｓ１１２のＮｏ）、処理をＳ１１３に進める。

Ｓ１１３において、更新部１４は、図７に示す第１例と同様の処理を実行する。Ｓ１１３に続いて、Ｓ１１４において、更新部１４は、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））の制約処理をする。具体的には、更新部１４は、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））を第３値に変更する。更新部１４は、Ｓ１１４を終えると、処理をＳ１１５に進める。

Ｓ１１５において、更新部１４は、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に、予め定められた演算により算出された加速値を加算する。例えば、加速値は、制約処理をする前の直前時刻における第２変数（ｙ´_ｉ（ｔ））に、所定係数（ｒ）と単位時間(Δｔ)との積（ｒΔｔ）を乗じた値である。更新部１４は、制約処理をする前の直前時刻における第２変数（ｙ´_ｉ（ｔ））を、Ｓ１０４において予め記憶している。例えば、更新部１４は、式（１９）に示す演算を実行する。

更新部１４は、Ｓ１１５を終えると、処理をＳ１１６に進める。更新部１４は、以上のようなＳ１１０とＳ１１６との間のループ処理をＮ回実行することにより、次の処理を実行する。

すなわち、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に基づき更新する。続いて、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値より小さい場合、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を第１値に変更する。また、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第２値より大きい場合、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を第２値に変更する。また、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が、第１値より小さいまたは第２値より大きい場合、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））を第３値に変更する。さらに、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に、制約処理をする前の対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に所定係数（ｒ）と単位時間（Δｔ）との積（ｒΔｔ）を乗じた加速値を加算する。

更新部１４は、Ｓ１１０とＳ１１６との間のループ処理をＮ回実行した場合、処理をＳ１０５に進める。

続いて、更新部１４は、Ｓ１０５とＳ１０９との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。

Ｓ１０６において、更新部１４は、Ｎ個の要素のそれぞれの対象時刻（ｔ＋Δｔ）における第１変数（ｘ_１（ｔ＋Δｔ）～ｘ_Ｎ（ｔ＋Δｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ＋Δｔ））を算出する。

ＱＵＢＯ問題の場合、更新部１４は、式（２０）を算出する。

また、ＨＯＢＯ問題の場合、更新部１４は、式（２１）を算出する。

続いて、Ｓ１０７において、更新部１４は、更新値（ｚ_ｉ（ｔ＋Δｔ））に係数（ｃ）を乗算することにより、外力（ｆ_ｉ（ｔ＋Δｔ））を算出する。具体的には、更新部１４は、式（２２）を算出する。

続いて、Ｓ１０８において、更新部１４は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δ））を、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に、外力（ｆ_ｉ（ｔ））と対象要素の対象時刻（ｔ＋Δ）における第１変数（ｘ_ｉ（ｔ＋Δ））とに基づく値に単位時間（Δｔ）を乗算した値を、加算することにより、算出する。具体的には、更新部１４は、式（２３）を算出する。

更新部１４は、以上のようなＳ１０５とＳ１０９との間のループ処理をＮ回実行することにより、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_１（ｔ＋Δｔ）～ｘ_Ｎ（ｔ＋Δｔ））に基づき更新する。

更新部１４は、Ｓ１０５とＳ１０９との間のループ処理をＮ回実行した場合、処理をＳ１１７に進める。

Ｓ１１７において、更新部１４は、直前時刻（ｔ）に単位時間（Δｔ）を加算して、対象時刻（ｔ＋Δｔ）を更新する。Ｓ１１８において、更新部１４は、Ｓ１０４からＳ１１７までの処理を、ｔが終了時刻（Ｔ）を超えるまで繰り返す。そして、更新部１４は、ｔが終了時刻（Ｔ）より大きくなった場合、本フローを終了する。

そして、更新部１４は、Ｎ個の要素のそれぞれについて、終了時刻（ｔ＝Ｔ）における第１変数（ｘ_ｉ（Ｔ））の符号に応じて、対応するスピンの値を算出する。

図８に示すフローチャートに従った処理を実行することにより、更新部１４は、Ｎ個の要素のそれぞれについて、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間毎に順次に、第１変数（ｘ_ｉ（ｔ＋Δｔ））および第２変数（ｙ_ｉ（ｔ＋Δｔ））を交互に更新する。また、図８に示すフローチャートに従った処理を実行することにより、更新部１４は、単位時間毎の更新処理において、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を算出した後に、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を算出する。さらに、図８に示すフローチャートに従った処理を実行することにより、更新部１４は、単位時間毎の更新処理において、Ｎ個の要素のそれぞれについて、第１変数（ｘ_ｉ（ｔ＋Δｔ））を更新した後、且つ、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を更新する前において、直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に加速値を加算する。

以上により、更新部１４は、シンプレクティック・オイラー法を用いて改良アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

なお、更新部１４は、Ｓ１１０～Ｓ１１６の処理を、Ｎ個の要素のそれぞれに対して実行する。更新部１４は、これらの処理をＭ_１個（Ｍ_１は、２以上、Ｎ以下の整数）の処理回路、プロセッサ、コンピュータまたはサーバ等の演算処理装置により並列に実行させてもよい。例えば、Ｍ_１は、Ｎの約数である。この場合、Ｎ個の要素のそれぞれは、Ｍ_１個の演算処理装置の何れか１つに対応している。そして、Ｍ_１個の演算処理装置のそれぞれは、Ｎ個の要素のうちの対応する要素について、Ｓ１１０～Ｓ１１６の処理を実行する。このような処理を実行することにより、更新部１４は、これらの処理を複数のＦＰＧＡ等の回路により分散して実行させて、高速に処理結果を得ることができる。

さらに、更新部１４は、Ｓ１０５～Ｓ１０９の処理も、Ｎ個の要素のそれぞれに対して実行する。更新部１４は、これらの処理をＭ_２個の処理回路、プロセッサ、コンピュータまたはサーバ等の演算処理装置により並列に実行させてもよい。例えば、Ｍ_２は、Ｎの約数である。この場合、Ｎ個の要素のそれぞれは、Ｍ_２個の演算処理装置の何れか１つに対応している。そして、Ｍ_２個の演算処理装置のそれぞれは、Ｎ個の要素のうちの対応する要素について、Ｓ１０５～Ｓ１０９の処理を実行する。このような処理を実行することにより、更新部１４は、これらの処理を複数のＦＰＧＡ等の回路により分散して実行させて、高速に処理結果を得ることができる。

（演算回路）
図９は、更新部１４を回路化する場合のブロック構成を示す図である。更新部１４は、回路化する場合、一例として図９に示すように構成される。なお、図９の説明において、直前時刻は、ｔ_１と表される。また、対象時刻は、ｔ_２と表される。

図９に示す更新部１４は、式（９）および式（１０）で表される改良アルゴリズムを実行する。更新部１４は、Ｘメモリ６６と、Ｙメモリ６７と、作用演算回路６８と、更新回路６９と、制御回路７０とを備える。

Ｘメモリ６６は、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｉ（ｔ_１））を記憶する。Ｘメモリ６６は、時刻の更新に伴って、記憶している直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｉ（ｔ_１））が書き換えられる。すなわち、Ｘメモリ６６は、対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｉ（ｔ_２））が算出された場合、算出された対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｉ（ｔ_２））が、新たな直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｉ（ｔ_１））として書き込まれる。制御回路７０は、演算に先だって、初期時刻におけるＮ個の第１変数ｘ_ｉをＸメモリ６６に書き込む。

Ｙメモリ６７は、直前時刻（ｔ_１）におけるＮ個の第２変数（ｙ_ｉ（ｔ_１））を記憶する。Ｙメモリ６７は、時刻の更新に伴って、記憶している直前時刻（ｔ_１）におけるＮ個の第２変数（ｙ_ｉ（ｔ_１））が書き換えられる。すなわち、Ｙメモリ６７は、対象時刻（ｔ_２）におけるＮ個の第２変数（ｙ_ｉ（ｔ_２））が算出された場合、算出された対象時刻（ｔ_２）におけるＮ個の第２変数（ｙ_ｉ（ｔ_２））が、新たな直前時刻（ｔ_１）におけるＮ個の第２変数（ｙ_ｉ（ｔ_１））として書き込まれる。制御回路７０は、演算に先だって、初期時刻におけるＮ個の第２変数ｙ_ｉをＹメモリ６７に書き込む。

作用演算回路６８は、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））をＸメモリ６６から取得する。そして、作用演算回路６８は、Ｎ個の要素のそれぞれについて、更新値（ｚ_ｉ（ｔ_１））を算出する。

更新回路６９は、Ｎ個の要素のそれぞれについて、更新値（ｚ_ｉ（ｔ_１））を作用演算回路６８から取得する。さらに、更新回路６９は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））をＸメモリ６６から取得し、直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））をＹメモリ６７から取得する。そして、更新回路６９は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出し、Ｘメモリ６６に記憶されている直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））を書き換える。これとともに、更新回路６９は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出し、Ｙメモリ６７に記憶されている直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））を書き換える。

制御回路７０は、対象時刻（ｔ_２）を単位時間（Δｔ）毎に順次に更新することにより、作用演算回路６８および更新回路６９に対して単位時間（Δｔ）毎の第１変数（ｘ_ｉ（ｔ））および第２変数（ｙ_ｉ（ｔ））を順次に算出させる。

さらに、制御回路７０は、１からＮまでのインクリメントすることによりインデックス（ｉ）を発生して、作用演算回路６８および更新回路６９に対して、Ｎ個の要素のそれぞれに対応する、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））および対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））をインデックス順に算出させる。なお、作用演算回路６８および更新回路６９は、複数個のインデックスに対応する複数個の第１変数（ｘ_ｉ（ｔ_２））および複数個の第２変数（ｙ_ｉ（ｔ_２））を並列に算出してもよい。

作用演算回路６８は、Ｊメモリ７１と、Ｈメモリ７２と、行列演算回路７３と、α関数回路７４と、第１加算回路７５とを有する。

Ｊメモリ７１は、（Ｎ×Ｎ）個の結合係数を含むＮ×Ｎの行列を記憶する。Ｊ_ｉ，_ｊは、行列に含まれるｉ行ｊ列の結合係数を表す。Ｊ_ｉ，_ｊは、組合せ最適化問題に対応するイジングモデルにおける、ｉ番目のスピンとｊ番目のスピンとの結合係数を表す。制御回路７０は、演算に先だって、ユーザ等により予め生成された行列をＪメモリ７１に書き込む。

Ｈメモリ７２は、Ｎ個の局所磁場係数を含む配列を記憶する。ｈ_ｉは、配列に含まれるｉ番目の局所磁場係数を表す。ｈ_ｉは、組合せ最適化問題に対応するイジングモデルにおける、ｉ番目のスピンに作用する局所磁場を表す。制御回路７０は、演算に先だって、ユーザ等により予め生成された配列をＨメモリ７２に書き込む。

行列演算回路７３は、Ｘメモリ６６から、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））を取得する。行列演算回路７３は、Ｎ個の要素のそれぞれについて、Ｊメモリ７１から対象の行に含まれるＮ個の結合係数Ｊ_ｉ，_ｊを取得する。そして、行列演算回路７３は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））と、対象の行に含まれるＮ個の結合係数Ｊ_ｉ，_ｊとの積和演算を実行する。

α関数回路７４は、Ｎ個の要素のそれぞれについて、Ｈメモリ７２から対象の局所磁場係数ｈ_ｉを取得する。α関数回路７４は、Ｎ個の要素のそれぞれについて、｛－ｈ_ｉα（ｔ_１）｝の演算を実行する。α（ｔ）は、予め設定された関数である。

第１加算回路７５は、Ｎ個の要素のそれぞれについて、行列演算回路７３による積和演算結果と、α関数回路７４による演算結果とを加算する。この演算により、第１加算回路７５は、Ｎ個の要素のそれぞれについて、式（２４）で表される直前時刻（ｔ_１）における更新値（ｚ_ｉ（ｔ_１））を出力する。

更新回路６９は、第１乗算回路７９と、Ｐ関数回路８０と、第２乗算回路８１と、第２加算回路８２と、第３乗算回路８３と、第３加算回路８４と、制約前Ｙメモリ８５と、第４乗算回路８６と、第４加算回路８７と、判定回路９１と、Ｘ制約回路９２と、Ｙ制約回路９３と、加速値算出回路９４と、加速値加算回路９５とを有する。

第１乗算回路７９は、Ｎ個の要素のそれぞれについて、更新値（ｚ_ｉ（ｔ_１））に係数である－ｃを乗算する。Ｐ関数回路８０は、Ｎ個の要素のそれぞれについて、｛－Ｄ＋ｐ（ｔ_１）｝の演算を実行する。第２乗算回路８１は、Ｎ個の要素のそれぞれについて、Ｘメモリ６６から直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））を取得する。そして、第２乗算回路８１は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））と、Ｐ関数回路８０の演算結果とを乗算する。

第２加算回路８２は、Ｎ個の要素のそれぞれについて、第１乗算回路７９の演算結果と、第２乗算回路８１の演算結果とを加算する。第３乗算回路８３は、第２加算回路８２の演算結果に単位時間であるΔｔを乗算する。

第３加算回路８４は、Ｎ個の要素のそれぞれについて、Ｙメモリ６７から直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））を取得する。第３加算回路８４は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））と、第３乗算回路８３の演算結果とを加算する。この演算により、第３加算回路８４は、Ｎ個の要素のそれぞれについて、式（２５）で表される対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を出力する。

そして、第３加算回路８４は、算出したＮ個の要素のそれぞれの対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を、制約前Ｙメモリ８５に書き込む。制約前Ｙメモリ８５は、Ｙ制約回路９３による制約前の、対象時刻（ｔ_２）におけるＮ個の第２変数（ｙ_ｉ（ｔ_２））を記憶する。

第４乗算回路８６は、Ｎ個の要素のそれぞれについて、制約前Ｙメモリ８５から対象時刻（ｔ_２）における、制約前の第２変数（ｙ_ｉ（ｔ_２））を取得する。第４乗算回路８６は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における制約前の第２変数（ｙ_ｉ（ｔ_２））に、｛ＤΔｔ｝を乗算する。

第４加算回路８７は、Ｎ個の要素のそれぞれについて、Ｘメモリ６６から直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））を取得する。第４加算回路８７は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））と、第４乗算回路８６の演算結果とを加算する。この演算により、第４加算回路８７は、Ｎ個の要素のそれぞれについて、式（２６）で表される対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を出力する。

判定回路９１は、Ｎ個の要素のそれぞれについて、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））が、予め定められた第１値以上、予め定められた第２値以下であるか否かを判断する。例えば、第１値は、－１である。第２値は、＋１である。判定回路９１は、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））が第１値より小さいまたは第２値より大きい場合、Ｘ制約回路９２およびＹ制約回路９３にイネーブル信号（ＥＮ）を与える。

Ｘ制約回路９２は、Ｎ個の要素のそれぞれについて、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を受け取る。Ｘ制約回路９２は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取らなかった場合、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を、そのままＸメモリ６６に書き込む。

Ｘ制約回路９２は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取った場合、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））に対して制約処理を実行し、制約処理をした後の第１変数（ｘ_ｉ（ｔ_２））をＸメモリ６６に書き込む。

Ｘ制約回路９２は、制約処理として、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））が第１値（例えば－１）より小さい場合、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を第１値（例えば－１）に変更する。Ｘ制約回路９２は、制約処理として、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））が第２値（例えば＋１）より大きい場合、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を第２値（例えば＋１）に変更する。

Ｙ制約回路９３は、Ｎ個の要素のそれぞれについて、制約前Ｙメモリ８５から、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を取得する。Ｙ制約回路９３は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取らなかった場合、制約前Ｙメモリ８５から取得した対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を、そのまま出力する。

Ｙ制約回路９３は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取った場合、制約前Ｙメモリ８５から取得した対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））に対して制約処理を実行し、制約処理をした後の第２変数（ｙ_ｉ（ｔ_２））を出力する。

ここで、Ｙ制約回路９３は、制約処理として、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を、第３値（例えば０）に変更する。Ｙ制約回路９３は、制約処理として、例えば、第１値以上、第２値以下のランダムな値に変更してもよい。

加速値算出回路９４は、Ｎ個の要素のそれぞれについて、Ｙメモリ６７から直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））を取得する。また、加速値算出回路９４は、所定係数（ｒ）を制御回路７０から演算に先だって受け取る。例えば、所定係数（ｒ）は、０より大きく１より小さい値である。そして、加速値算出回路９４は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））と、所定係数（ｒ）と単位時間（Δｔ）との積（ｒΔｔ）とを乗算して加速値を算出する。

加速値加算回路９５は、Ｎ個の要素のそれぞれについて、Ｙ制約回路９３から、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を取得する。また、加速値加算回路９５は、Ｎ個の要素のそれぞれについて、加速値算出回路９４から加速値（ｒΔｔ×（ｙ_ｉ（ｔ_１）））を取得する。加速値加算回路９５は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））に加速値（ｒΔｔ×（ｙ_ｉ（ｔ_１）））を加算する。そして、加速値加算回路９５は、加速値（ｒΔｔ×（ｙ_ｉ（ｔ_１）））が加算された対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を、Ｙメモリ６７に書き込む。

以上により、更新部１４は、式（９）の改良アルゴリズムを実行して、終了時刻（Ｔ）におけるＮ個の第１変数（ｘ_ｉ（Ｔ））およびＮ個の第２変数（ｙ_ｉ（Ｔ））を算出することができる。

なお、上述した更新部１４は、単位時間毎に、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出した後に、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出する。これに代えて、更新部１４は、単位時間毎に、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出した後に、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出してもよい。

この場合、制約前Ｙメモリ８５は、Ｙ制約回路９３による制約前の、直前時刻（ｔ_１）におけるＮ個の第２変数（ｙ_ｉ（ｔ_１））を記憶する。そして、第４加算回路８７は、Ｎ個の要素のそれぞれについて、式（２７）で表される対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を出力する。

行列演算回路７３は、Ｘメモリ６６から、対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｊ（ｔ_２））を取得する。そして、第１加算回路７５は、Ｎ個の要素のそれぞれについて、式（２８）で表される対象時刻（ｔ_２）における更新値（ｚ_ｉ（ｔ_２））を出力する。

また、第１乗算回路７９は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における更新値（ｚ_ｉ（ｔ_２））に－ｃを乗算する。Ｐ関数回路８０は、Ｎ個の要素のそれぞれについて、｛－Ｄ＋ｐ（ｔ_２）｝の演算を実行する。第２乗算回路８１は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））と、Ｐ関数回路８０の演算結果とを乗算する。

そして、第３加算回路８４は、Ｎ個の要素のそれぞれについて、式（２９）で表される対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を出力する。

このような処理によっても、更新部１４は、改良アルゴリズムを実行して、終了時刻（Ｔ）におけるＮ個の第１変数（ｘ_ｉ（Ｔ））およびＮ個の第２変数（ｙ_ｉ（Ｔ））を算出することができる。

なお、更新回路６９は、Ｎ個の要素のそれぞれに対して、第２変数（ｙ_ｉ（ｔ_１））に基づき第１変数（ｘ_ｉ（ｔ_１））を更新するＸ更新処理、第１変数（ｘ_ｉ（ｔ_２））を第１値または第２値に変更するＸ制約処理、第２変数（ｙ_ｉ（ｔ_２））を第３値に変更するＹ制約処理、および、第２変数（ｙ_ｉ（ｔ_２））に対して加速値（ｒΔｔ×（ｙ_ｉ（ｔ_１）））を加算する加速値加算処理を実行する。更新回路６９は、これらの処理をＭ_１個の処理回路により並列に実行してもよい。この場合、Ｎ個の要素のそれぞれは、Ｍ_１個の処理回路の何れか１つに対応している。そして、Ｍ_１個の処理回路のそれぞれは、Ｎ個の要素のうちの対応する要素について、Ｘ更新処理、Ｘ制約処理、Ｙ制約処理および加速値加算処理を実行する。このような処理を実行することにより、更新回路６９は、これらの処理を複数の処理回路により分散して実行させて、高速に処理結果を得ることができる。

さらに、作用演算回路６８は、Ｎ個の要素のそれぞれについて、更新値（ｚ_ｉ（ｔ_１））を算出する。作用演算回路６８は、更新値（ｚ_ｉ（ｔ_１））の算出処理を、Ｍ_２個の処理回路により並列に実行させてもよい。この場合、Ｎ個の要素のそれぞれは、Ｍ_２個の処理回路の何れか１つに対応している。そして、Ｍ_２個の処理回路のそれぞれは、Ｎ個の要素のうちの対応する要素について、更新値（ｚ_ｉ（ｔ_１））の算出処理を実行する。また、作用演算回路６８を並列化する場合、行列演算回路７３の並列度と、α関数回路７４の並列度と、第１加算回路７５の並列度とが異なっていてもよいし、行列演算回路７３のみが並列化されていてもよい。このような処理を実行することにより、作用演算回路６８は、更新値（ｚ_ｉ（ｔ_１））の算出処理を複数の処理回路により分散して実行させて、高速に処理結果を得ることができる。特に、行列演算回路７３が並列化されることにより、作用演算回路６８は、行列演算に伴う膨大な量の乗算および加算処理を効率良く実行することができる。

（システム構成）
図１０は、情報処理システム１００の構成を示す図である。改良アルゴリズムは、例えば図１０に示すような情報処理システム１００により実行させることが可能である。これにより、情報処理システム１００は、大規模な組合せ最適化問題を、並列処理により高速に解くことができる。

情報処理システム１００は、管理サーバ１０１と、ネットワーク１０２と、複数の計算サーバ１０３（１０３ａ～１０３ｃ）（情報処理装置）と、複数のケーブル１０４（１０４ａ～１０４ｃ）と、スイッチ１０５を備える。また、端末装置１０６は、情報処理システム１００と通信可能である。管理サーバ１０１、複数の計算サーバ１０３（１０３ａ～１０３ｃ）、端末装置１０６は、ネットワーク１０２を介して互いにデータ通信をする。ネットワーク１０２は、例えば、複数のコンピュータネットワークが相互に接続されたインターネットである。ネットワーク１０２は、通信媒体として有線、無線、または、これらの組合せであってよい。

また、複数の計算サーバ１０３（１０３ａ～１０３ｃ）は、それぞれケーブル１０４（１０４ａ～１０４ｃ）を介してスイッチ１０５に接続される。複数のケーブル１０４（１０４ａ～１０４ｃ）およびスイッチ１０５は、計算サーバ１０３と計算サーバ１０３と間のインターコネクトを形成する。複数の計算サーバ１０３（１０３ａ～１０３ｃ）のそれぞれは、インターコネクトを介して互いにデータ通信をすることも可能である。スイッチ１０５は、例えば、Ｉｎｆｉｎｉｂａｎｄのスイッチでる。ケーブル１０４ａ～１０４ｃは、例えば、Ｉｎｆｉｎｉｂａｎｄのケーブルである。ただし、スイッチ１０５およびケーブル１０４は、有線ＬＡＮのスイッチ／ケーブルであってもよい。ケーブル１０４およびスイッチ１０５で使われる通信規格および通信プロトコルについては、特に問わない。端末装置１０６は、例えば、ノートＰＣ、デスクトップＰＣ、スマートフォン、タブレットまたは車載端末装置である。

改良アルゴリズムを用いた組合せ最適化問題の求解処理は、並列的および分散的に実行可能である。従って、複数の計算サーバ１０３（１０３ａ～１０３ｃ）のそれぞれおよび／または計算サーバ１０３（１０３ａ～１０３ｃ）のプロセッサは、改良アルゴリズムを用いた組合せ最適化問題の求解処理における一部の計算処理の一部のステップを分担して実行してもよいし、異なる変数に対する同一の計算処理を並列的に実行してもよい。この場合、管理サーバ１０１は、例えば、ユーザによって入力された組合せ最適化問題を各計算サーバ１０３に処理可能な形式に変換し、各計算サーバ１０３に実行させる。そして、管理サーバ１０１は、各計算サーバ１０３から計算結果を取得し、集約した計算結果を組合せ最適化問題の解に変換する。

図１０には、３台の計算サーバ１０３（１０３ａ～１０３ｃ）が示されている。しかし、情報処理システム１００に含まれる計算サーバ１０３の台数は、３台に限定されない。例えば、情報処理システム１００は、１台であっても、２台であっても、４台以上であってもよい。また、情報処理システム１００は、含んでいる複数の計算サーバ１０３のうちの一部を用いて、組合せ最適化問題の求解の処理を実行してもよい。計算サーバ１０３は、どのような種類の情報処理装置であってもよい。例えば、計算サーバ１０３は、データセンターに設置されたサーバであってもよいし、オフィスに設置されたデスクトップＰＣであってもよい。また、計算サーバ１０３は、異なるローケーションに設置された複数の種類のコンピュータであってもよい。例えば、計算サーバ１０３は、汎用的なコンピュータであってもよいし、専用の電子回路または、これらの組合せであってもよい。

図１１は、管理サーバ１０１の構成を示す図である。管理サーバ１０１は、例えば、中央演算処理装置（ＣＰＵ）とメモリとを含むコンピュータである。管理サーバ１０１は、プロセッサ１１０と、記憶部１１４と、通信回路１１５と、入力回路１１６と、出力回路１１７とを備える。プロセッサ１１０、記憶部１１４、通信回路１１５、入力回路１１６、出力回路１１７は、互いにバス１２０を介して接続される。プロセッサ１１０は、内部の機能構成として、管理部１１１と、変換部１１２と、制御部１１３を含む。

プロセッサ１１０は、演算を実行し、管理サーバ１０１の制御を行う電子回路である。プロセッサ１１０は、例えば、ＣＰＵ、マイクロプロセッサ、ＡＳＩＣ、ＦＰＧＡ、ＰＬＤまたはこれらの組合せであってよい。管理部１１１は、ユーザの端末装置１０６を介して管理サーバ１０１の操作を行うためのインタフェースを提供する。管理部１１１が提供するインタフェースは、例えば、ＡＰＩ、ＣＬＩまたはウェブページ等である。例えば、ユーザは、管理部１１１を介して組合せ最適化問題の情報の入力を行ったり、計算された組合せ最適化問題の解の閲覧および／またはダウンロードを行ったりする。変換部１１２は、組合せ最適化問題に関するパラメータを入力して、入力したパラメータを各計算サーバ１０３が処理可能な形式に変換する。制御部１１３は、各計算サーバ１０３に制御指令を送信する。制御部１１３が各計算サーバ１０３から計算結果を取得した後、変換部１１２は、複数の計算結果を集約し、組合せ最適化問題の解に変換し、組合せ最適化問題の解を出力する。

記憶部１１４は、管理サーバ１０１のプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。プログラムは、ＯＳとアプリケーションの両方を含む。記憶部１１４は、揮発性メモリ、不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリは、例えば、ＤＲＡＭまたはＳＲＡＭ等である。不揮発性メモリは、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭまたはＭＲＡＭ等である。また、記憶部１１４は、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置等であってもよい。

通信回路１１５は、ネットワーク１０２に接続された各装置との間でデータの送受信を行う。通信回路１１５は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ Ｉｎｔｅｒｆａｃｅ Ｃａｒｄ）である。ただし、通信回路１１５は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。入力回路１１６は、管理サーバ１０１へのデータ入力を実現する。入力回路１１６は、外部ポートとして、例えば、ＵＳＢ、ＰＣＩ－Ｅｘｐｒｅｓｓなどを含む。操作装置１１８は、入力回路１１６に接続される。操作装置１１８は、ユーザが管理サーバ１０１に情報を入力するための装置である。操作装置１１８は、例えば、キーボード、マウス、タッチパネル、音声認識装置などであるが、これらに限られない。出力回路１１７は、管理サーバ１０１からのデータ出力を実現する。出力回路１１７は、外部ポートとしてＨＤＭＩ（登録商標）またはＤｉｓｐｌａｙＰｏｒｔ等を含む。例えば、表示装置１１９は、出力回路１１７に接続される。表示装置１１９は、例えば、ＬＣＤ（液晶ディスプレイ）、有機ＥＬ（有機エレクトロルミネッセンス）ディスプレイまたはプロジェクタ等であるが、これらに限られない。

管理者は、操作装置１１８および表示装置１１９を使って、管理サーバ１０１のメンテナンスを行うことができる。なお、操作装置１１８および表示装置１１９は、管理サーバ１０１に組み込まれたものであってもよい。また、操作装置１１８および表示装置１１９は、管理サーバ１０１に接続されていなくてもよい。例えば、管理者は、ネットワーク１０２と通信可能な端末装置を用いて管理サーバ１０１のメンテナンスを行ってもよい。

図１２は、管理サーバ１０１の記憶部１１４に記憶されるデータを示す図である。記憶部１１４は、例えば、問題データ１１４Ａと、計算データ１１４Ｂと、管理プログラム１１４Ｃと、変換プログラム１１４Ｄと、制御プログラム１１４Ｅとを記憶する。例えば、問題データ１１４Ａは、組合せ最適化問題のデータを含む。例えば、計算データ１１４Ｂは、各計算サーバ１０３から収集された計算結果を含む。例えば、管理プログラム１１４Ｃは、管理部１１１の機能を実現するプログラムである。例えば、変換プログラム１１４Ｄは、変換部１１２の機能を実現するプログラムである。例えば、制御プログラム１１４Ｅは、制御部１１３の機能を実現するプログラムである。

図１３は、計算サーバ１０３ａの構成を示す図である。他の計算サーバ１０３は、計算サーバ１０３ａと同様の構成であってもよいし、計算サーバ１０３ａと異なる構成であってもよい。

計算サーバ１０３ａは、例えば、通信回路１３１と、共有メモリ１３２と、プロセッサ１３３ａ～１３３ｄと、ストレージ１３４と、ホストバスアダプタ１３５とを備える。通信回路１３１、共有メモリ１３２、プロセッサ１３３ａ～１３３ｄ、ストレージ１３４、ホストバスアダプタ１３５は、バス１３６を介して互いに接続される。

通信回路１３１は、ネットワーク１０２に接続された各装置との間でデータの送受信を行う。通信回路１３１は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ Ｉｎｔｅｒｆａｃｅ Ｃａｒｄ）である。ただし、通信回路１３１は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。共有メモリ１３２は、プロセッサ１３３ａ～１３３ｄからアクセス可能なメモリである。共有メモリ１３２は、例えば、ＤＲＡＭおよびＳＲＡＭ等の揮発性メモリである。共有メモリ１３２は、不揮発性メモリ等の他の種類のメモリを含んでもよい。プロセッサ１３３ａ～１３３ｄは、共有メモリ１３２を介してデータを共有する。なお、共有メモリ１３２は、計算サーバ１０３ａの全てメモリにより構成されていなくてもよい。例えば、計算サーバ１０３ａの一部のメモリは、プロセッサ１３３ａ～１３３ｄのうちのいずれかからのみからアクセスできるローカルメモリであってもよい。

プロセッサ１３３ａ～１３３ｄは、計算処理を実行する電子回路である。プロセッサ１３３ａ～１３３ｄは、例えば、ＣＰＵ、ＧＰＵ、ＦＰＧＡ、ＡＳＩＣのいずれであってもよいし、これらの組合せであってもよい。また、プロセッサ１３３ａ～１３３ｄは、ＣＰＵコアまたはＣＰＵスレッドであってもよい。プロセッサ１３３ａ～１３３ｄがＣＰＵである場合、計算サーバ１０３ａが備えるソケット数については、特に問わない。また、プロセッサ１３３ａ～１３３ｄは、ＰＣＩ ｅｘｐｒｅｓｓなどのバスを介して計算サーバ１０３ａのその他の構成要素に接続されていてもよい。

図１３の例では、計算サーバ１０３ａは、４つのプロセッサ１３３ａ～１３３ｄを備える。しかし、１台の計算サーバ１０３ａに含まれるプロセッサ数は、４個に限られない。

ストレージ１３４は、計算サーバ１０３ａのプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、ＯＳとアプリケーションの両方を含むものとする。ストレージ１３４は、揮発性メモリおよび不揮発性メモリ、またはこれらの組合せであってもよい。揮発性メモリは、たとえば、ＤＲＡＭまたはＳＲＡＭである。不揮発性メモリは、例えば、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭまたはＭＲＡＭ等である。また、ストレージ１３４は、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置を含んでもよい。

ホストバスアダプタ１３５は、他の計算サーバ１０３との間のデータ通信を実現する。ホストバスアダプタ１３５は、ケーブル１０４ａを介してスイッチ１０５に接続されている。ホストバスアダプタ１３５は、例えば、ＨＣＡ（Ｈｏｓｔ Ｃｈａｎｎｅｌ Ａｄａｐｔｏｒ）である。ホストバスアダプタ１３５、ケーブル１０４ａ、スイッチ１０５で高スループットを実現可能なインターコネクトを形成することにより、並列的な計算処理の速度を向上させることができる。

図１４は、ストレージ１３４に記憶されるデータを示す図である。ストレージ１３４は、例えば、計算データ１３４Ａと、計算プログラム１３４Ｂと、制御プログラム１３４Ｃとを記憶する。計算データ１３４Ａは、計算サーバ１０３ａの計算途中のデータまたは計算結果を含む。計算データ１３４Ａの少なくとも一部は、共有メモリ１３２、プロセッサのキャッシュまたはプロセッサのレジスタ等の、異なる記憶階層に記憶されてもよい。計算プログラム１３４Ｂは、各プロセッサ１３３における計算処理および、共有メモリ１３２およびストレージ１３４へのデータの保存処理を実現するプログラムである。制御プログラム１３４Ｃは、管理サーバ１０１の制御部１１３から送信された指令に基づき、計算サーバ１０３ａを制御し、計算サーバ１０３ａの計算結果を管理サーバ１０１に送信するプログラムである。

このような計算サーバ１０３ａは、プロセッサ１３３ａ～１３３ｄが組合せ最適化問題を解くためのプログラムを実行する。このプログラムは、計算サーバ１０３ａを、入力部１２、更新部１４および出力部１６として機能させる。

このような構成の情報処理システム１００は、ＰＣクラスタとして利用することが可能である。ＰＣクラスタは、複数台のコンピュータを接続し、１台のコンピュータでは得られない計算性能を実現するシステムである。情報処理システム１００は、例えばＭＰＩ（Ｍｅｓｓａｇｅ Ｐａｓｓｉｎｇ Ｉｎｔｅｒｆａｃｅ）を使うことにより、複数の計算サーバ１０３にメモリが分散して配置されている構成であっても、並列的な計算を実行することが可能である。

なお、情報処理システム１００は、高速リンクで接続された複数のＧＰＵであってもよい。この場合、複数のＧＰＵのそれぞれは、計算サーバ１０３と同様の処理を実行する。

本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組合せにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組合せてもよい。

（技術案）
なお、上記の実施形態を、以下の技術案にまとめることができる。

技術案１
組合せ最適化問題を解く計算装置であって、
第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力部と、
を備え、
前記複数の要素は、前記組合せ最適化問題の複数の離散変数に対応し、
前記第１変数および前記第２変数のそれぞれは、実数により表され、
前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第１変数を前記第２変数に基づき更新し、
前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値に変更し、前記第２変数を予め定められた第３値に変更し、
前記第１変数が、前記第１値より大きい予め定められた第２値より大きい場合、前記第１変数を前記第２値に変更し、前記第２変数を前記第３値に変更し、
前記第２変数に、予め定められた演算により算出された加速値を加算する
計算装置。

技術案２
前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第１変数を更新する前において、前記第２変数を更新し、
前記第１変数を更新した後において、前記第２変数に前記加速値を加算する
技術案１に記載の計算装置。

技術案３
前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第１変数を更新した後において、前記第２変数を更新し、
前記第１変数を更新した後、且つ、前記第２変数を更新する前において、前記第２変数に前記加速値を加算する
技術案１に記載の計算装置。

技術案４
前記出力部は、
前記終了時刻における前記複数の要素のそれぞれについて、前記第１変数を予め設定されたしきい値により２値化した離散変数の値を算出し、
算出した前記複数の離散変数の値を前記組合せ最適化問題の解として出力する
技術案１から３の何れか１項に記載の計算装置。

技術案５
前記加速値は、直前時刻における対応する前記第２変数に、所定係数と前記単位時間との積を乗じた値である
技術案１から４の何れか１項に記載の計算装置。

技術案６
前記所定係数は、予め定められた定数である
技術案５に記載の計算装置。

技術案７
前記所定係数は、乱数に基づき定まる値である
技術案５に記載の計算装置。

技術案８
前記第１値は、－１、
前記第２値は、＋１、および、
前記第３値は、０である
技術案１から７の何れか１項に記載の計算装置。

技術案９
前記複数の要素のそれぞれは、複数の処理回路の何れか１つに対応しており、
前記複数の処理回路のそれぞれは、前記複数の要素のうちの対応する要素について、前記第１変数を前記第２変数に基づき更新する処理、前記第１変数を前記第１値または前記第２値に変更する処理および前記第２変数を前記第３値に変更する処理、および、前記第２変数に前記加速値を加算する処理を実行する
技術案１から８何れか１項に記載の計算装置。

技術案１０
前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、対象時刻における前記第１変数を、前記対象時刻より単位時間前の直前時刻における前記第１変数に、前記第２変数と予め定められた定数と前記単位時間とを乗算した値を加算することにより、算出する
技術案１から９の何れか１項に記載の計算装置。

技術案１１
前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、
前記複数の要素のそれぞれの前記第１変数と、対象要素と前記複数の要素のそれぞれとの組毎に予め定められる作用係数とに基づき、外力を算出し、
前記直前時刻における前記第２変数に、前記外力と前記第１変数とにより定まる値に前記単位時間とを乗算した値を加算することにより、前記対象時刻における前記第２変数を算出する
技術案１０に記載の計算装置。

技術案１２
前記組合せ最適化問題は、Ｎ個の離散変数を含み、
前記更新部は、前記Ｎ個の離散変数のうちのｉ番目の離散変数に対応するｉ番目の要素の、前記対象時刻における前記第１変数を式（１０１）または（１０２）により算出し、

Ｎは、２以上の整数であり、
ｉは、１からＮまでの任意の整数であり、
Ｄは、前記定数であり、
Δｔは、前記単位時間であり、
ｔは、前記直前時刻であり、
ｔ＋Δｔは、前記対象時刻であり、
ｘ_ｉ（ｔ）は、前記ｉ番目の要素の前記直前時刻における前記第１変数であり、
ｙ_ｉ（ｔ）は、前記ｉ番目の要素の前記直前時刻における前記第２変数であり、
ｘ_ｉ（ｔ＋Δｔ）は、前記ｉ番目の要素の前記対象時刻における前記第１変数であり、
ｙ_ｉ（ｔ＋Δｔ）は、前記ｉ番目の要素の前記対象時刻における前記第２変数である
技術案１１に記載の計算装置。

技術案１３
前記組合せ最適化問題は、ＱＵＢＯ（Ｑｕａｄｒａｔｉｃ Ｕｎｃｏｎｓｔｒａｉｎｅｄ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）問題であり、
前記更新部は、前記ｉ番目の要素の前記対象時刻における前記第２変数を式（１０３）または式（１０４）により算出し、

ｆ_ｉ（ｔ＋Δｔ）は、式（１０５）により表され、ｆ_ｉ（ｔ）は、式（１０６）により表され、

ｚ_ｉ（ｔ＋Δｔ）は、式（１０７）により表され、
ｚ_ｉ（ｔ）は、式（１０８）により表され、

ｊは、１からＮまでの任意の整数であり、
ｈ_ｉは、Ｎ個の局所磁場係数を含む予め定められた配列に含まれる、ｉ番目の局所磁場係数であり、
Ｊ_ｉ，ｊは、Ｎ×Ｎ個の結合係数を含む予め定められた行列に含まれる、ｉ行、ｊ列の結合係数であり、
ｃは、係数であり、
ｘ_ｊ（ｔ）は、前記Ｎ個の離散変数のうちのｊ番目の離散変数に対応するｊ番目の要素の前記直前時刻における前記第１変数であり、
ｘ_ｊ（ｔ＋Δｔ）は、前記ｊ番目の要素の前記対象時刻における前記第１変数であり、
ｐ（ｔ）は、ｔを変数とする予め定められた関数であり、ｔが増大するに従って増加し、ｔが前記初期時刻において０となり、ｔが前記終了時刻において１となり、
α（ｔ）は、ｔを変数とする予め定められた関数である
技術案１２に記載の計算装置。

技術案１４
情報処理装置により、組合せ最適化問題を解くための計算方法であって、
前記情報処理装置により、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新ステップと、
前記情報処理装置により、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力ステップと、
を実行し、
前記複数の要素は、前記組合せ最適化問題の複数の離散変数に対応し、
前記第１変数および前記第２変数のそれぞれは、実数により表され、
前記更新ステップにおける前記単位時間毎の更新処理において、前記情報処理装置は、前記複数の要素のそれぞれについて、
前記第１変数を前記第２変数に基づき更新し、
前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値に変更し、前記第２変数を予め定められた第３値に変更し、
前記第１変数が、前記第１値より大きい予め定められた第２値より大きい場合、前記第１変数を前記第２値に変更し、前記第２変数を前記第３値に変更し、
前記第２変数に、予め定められた演算により算出された加速値を加算する
計算方法。

技術案１５
情報処理装置を、組合せ最適化問題を解く計算装置として機能させるためのプログラムであって、
前記情報処理装置を、
第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力部と、
して機能させ、
前記複数の要素は、前記組合せ最適化問題の複数の離散変数に対応し、
前記第１変数および前記第２変数のそれぞれは、実数により表され、
前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第１変数を前記第２変数に基づき更新し、
前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値に変更し、前記第２変数を予め定められた第３値に変更し、
前記第１変数が、前記第１値より大きい予め定められた第２値より大きい場合、前記第１変数を前記第２値に変更し、前記第２変数を前記第３値に変更し、
前記第２変数に、予め定められた演算により算出された加速値を加算する
プログラム。

１０ 計算装置
１２ 入力部
１４ 更新部
１６ 出力部
１００ 情報処理システム
１０１ 管理サーバ
１０２ ネットワーク
１０３，１０３ａ，１０３ｂ，１０３ｃ 計算サーバ
１０４，１０４ａ，１０４ｂ，１０４ｃ ケーブル
１０５ スイッチ
１０６ 端末装置
１１０ プロセッサ
１１１ 管理部
１１２ 変換部
１１３ 制御部
１１４ 記憶部
１１５ 通信回路
１１６ 入力回路
１１７ 出力回路
１２０ バス
１３１ 通信回路
１３２ 共有メモリ
１３３ａ，１３３ｂ，１３３ｃ，１３３ｄ プロセッサ
１３４ ストレージ
１３５ ホストバスアダプタ

Claims (9)

1. 組合せ最適化問題を解く計算装置であって、
   第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
   前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力部と、
   を備え、
   前記複数の要素は、前記組合せ最適化問題の複数の離散変数に対応し、
   前記第１変数および前記第２変数のそれぞれは、実数により表され、
   前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
   前記第１変数を前記第２変数に基づき更新し、
   前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値に変更し、前記第２変数を予め定められた第３値に変更し、
   前記第１変数が、前記第１値より大きい予め定められた第２値より大きい場合、前記第１変数を前記第２値に変更し、前記第２変数を前記第３値に変更し、
   前記第２変数に、予め定められた演算により算出された加速値を加算し、
   前記加速値は、直前時刻における対応する前記第２変数に、所定係数と前記単位時間との積を乗じた値である
   計算装置。
2. 前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
   前記第１変数を更新する前において、前記第２変数を更新し、
   前記第１変数を更新した後において、前記第２変数に前記加速値を加算する
   請求項１に記載の計算装置。
3. 前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
   前記第１変数を更新した後において、前記第２変数を更新し、
   前記第１変数を更新した後、且つ、前記第２変数を更新する前において、前記第２変数に前記加速値を加算する
   請求項１に記載の計算装置。
4. 前記出力部は、
   前記終了時刻における前記複数の要素のそれぞれについて、前記第１変数を予め設定されたしきい値により２値化した離散変数の値を算出し、
   算出した前記複数の離散変数の値を前記組合せ最適化問題の解として出力する
   請求項１から３の何れか１項に記載の計算装置。
5. 前記複数の要素のそれぞれは、複数の処理回路の何れか１つに対応しており、
   前記複数の処理回路のそれぞれは、前記複数の要素のうちの対応する要素について、前記第１変数を前記第２変数に基づき更新する処理、前記第１変数を前記第１値または前記第２値に変更する処理および前記第２変数を前記第３値に変更する処理、および、前記第２変数に前記加速値を加算する処理を実行する
   請求項１から４の何れか１項に記載の計算装置。
6. 前記単位時間毎の更新処理において、前記更新部は、
   前記複数の要素のそれぞれについて、対象時刻における前記第１変数を、前記対象時刻より単位時間前の直前時刻における前記第１変数に、前記第２変数と予め定められた定数と前記単位時間とを乗算した値を加算することにより、算出する
   請求項１から５の何れか１項に記載の計算装置。
7. 前記単位時間毎の更新処理において、前記更新部は、
   前記複数の要素のそれぞれについて、
   前記複数の要素のそれぞれの前記第１変数と、対象要素と前記複数の要素のそれぞれとの組毎に予め定められる作用係数とに基づき、外力を算出し、
   前記直前時刻における前記第２変数に、前記外力と前記第１変数とにより定まる値に前記単位時間とを乗算した値を加算することにより、前記対象時刻における前記第２変数を算出する
   請求項６に記載の計算装置。
8. 情報処理装置により、組合せ最適化問題を解くための計算方法であって、
   前記情報処理装置により、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新ステップと、
   前記情報処理装置により、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力ステップと、
   を実行し、
   前記複数の要素は、前記組合せ最適化問題の複数の離散変数に対応し、
   前記第１変数および前記第２変数のそれぞれは、実数により表され、
   前記更新ステップにおける前記単位時間毎の更新処理において、前記情報処理装置は、前記複数の要素のそれぞれについて、
   前記第１変数を前記第２変数に基づき更新し、
   前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値に変更し、前記第２変数を予め定められた第３値に変更し、
   前記第１変数が、前記第１値より大きい予め定められた第２値より大きい場合、前記第１変数を前記第２値に変更し、前記第２変数を前記第３値に変更し、
   前記第２変数に、予め定められた演算により算出された加速値を加算し、
   前記加速値は、直前時刻における対応する前記第２変数に、所定係数と前記単位時間との積を乗じた値である
   計算方法。
9. 情報処理装置を、組合せ最適化問題を解く計算装置として機能させるためのプログラムであって、
   前記情報処理装置を、
   第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
   前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力部と、
   して機能させ、
   前記複数の要素は、前記組合せ最適化問題の複数の離散変数に対応し、
   前記第１変数および前記第２変数のそれぞれは、実数により表され、
   前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
   前記第１変数を前記第２変数に基づき更新し、
   前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値に変更し、前記第２変数を予め定められた第３値に変更し、
   前記第１変数が、前記第１値より大きい予め定められた第２値より大きい場合、前記第１変数を前記第２値に変更し、前記第２変数を前記第３値に変更し、
   前記第２変数に、予め定められた演算により算出された加速値を加算し、
   前記加速値は、直前時刻における対応する前記第２変数に、所定係数と前記単位時間との積を乗じた値である
   プログラム。

Images