ホーム > 特許ランキング > 日立ヴァンタラ株式会社
知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2025-02-03
(45)【発行日】2025-02-12
(54)【発明の名称】情報処理システム、情報処理方法、及び情報処理プログラム
(51)【国際特許分類】
G06F 17/10 20060101AFI20250204BHJP
G06N 99/00 20190101ALI20250204BHJP
【ＦＩ】
G06F17/10 Z
G06N99/00 180
【請求項の数】
14
(21)【出願番号】P 2021062592
(22)【出願日】2021-04-01
(65)【公開番号】P2022158010
(43)【公開日】2022-10-14
【審査請求日】2024-03-07
(73)【特許権者】
【識別番号】524132520
【氏名又は名称】日立ヴァンタラ株式会社
(74)【代理人】
【識別番号】110002365
【氏名又は名称】弁理士法人サンネクスト国際特許事務所
(72)【発明者】
【氏名】奥山 拓哉
【審査官】坂東 博司
(56)【参考文献】
【文献】国際公開第２０２０／２０２３１２（ＷＯ，Ａ１）
【文献】特開２０２０－２０４９２９（ＪＰ，Ａ）
【文献】米国特許出願公開第２０１６／０２６００１３（ＵＳ，Ａ１）
【文献】OKUYAMA, Takuya et al.,，"Binary optimization by momentum annealing"， Physical Review E [online],，Vol. 100，Iss. 1，米国，American Physical Society，2019年07月10日，pp.012111-1～012111-9，［online］，［令和 6年11月22日検索］，インターネット＜URL：https://journals.aps.org/pre/pdf/10.1103/PhysRevE.100.012111＞, DOI:10.1103/PhysRevE.100.012111，DOI:10.1103/PhysRevE.100.012111
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｆ １７／１０
Ｇ０６Ｎ ９９／００
(57)【特許請求の範囲】
【請求項1】
行列積を用いて最適化問題の最適解探索を行う情報処理システムであって、記憶部と協働して処理を実行する処理部を有し、
前記処理部は、
前記最適化問題を、第１の状態変数ベクトルと所定行列と第２の状態変数ベクトルとの積を含んだ目的関数で表される二次計画問題に変換する変換処理と、
前記所定行列と前記第１の状態変数ベクトルとを乗算して第１ベクトルを計算し、該第１ベクトルに基づくパラメータを含んだ確率分布に基づいて前記第２の状態変数ベクトルの確率的更新を行い、前記所定行列と前記第２の状態変数ベクトルとを乗算して第２ベクトルを計算し、該第２ベクトルに基づくパラメータを含んだ確率分布に基づいて前記第１の状態変数ベクトルの確率的更新を行う処理を繰り返すことで、前記第１の状態変数ベクトル及び前記第２の状態変数ベクトルの確率的更新を行う状態更新処理と、
前記第１の状態変数ベクトルと前記第２ベクトルとの内積によって前記積を計算することで、前記目的関数の関数値を算出する関数値算出処理と
を実行することを特徴とする情報処理システム。
【請求項2】
請求項１に記載の情報処理システムであって、前記処理部は、
前記状態更新処理及び前記関数値算出処理を、前記第１の状態変数ベクトル及び前記第２の状態変数ベクトルに対応する複数のレプリカ毎に独立かつ並行に実行する
ことを特徴とする情報処理システム。
【請求項3】
請求項２に記載の情報処理システムであって、複数の演算装置を有し、
前記処理部は、
前記レプリカ毎の前記状態更新処理及び前記関数値算出処理を、前記複数の演算装置によって実行する
ことを特徴とする情報処理システム。
【請求項4】
請求項３に記載の情報処理システムであって、前記処理部は、
前記関数値算出処理によって前記レプリカ毎の前記関数値を算出し、
前記レプリカ毎の前記関数値及び逆温度の差分に基づいて前記レプリカ毎の重みを算出し、該重みに基づいて、前記複数のレプリカから選択するレプリカの選択数を決定するリサンプリング処理を実行する
ことを特徴とする情報処理システム。
【請求項5】
請求項４に記載の情報処理システムであって、前記選択数は、前記複数のレプリカの合計数と同一である
ことを特徴とする情報処理システム。
【請求項6】
請求項４に記載の情報処理システムであって、前記選択数は、前記複数のレプリカの合計数よりも多い
ことを特徴とする情報処理システム。
【請求項7】
請求項４に記載の情報処理システムであって、前記選択数は、前記複数のレプリカの合計数よりも少ない
ことを特徴とする情報処理システム。
【請求項8】
請求項４に記載の情報処理システムであって、前記処理部は、
前記リサンプリング処理において、前記選択数に応じて、前記選択するレプリカが配置される前記演算装置に該レプリカのコピーを配置可能である場合には該演算装置を該レプリカのコピーを配置するコピー先として決定し、配置不可能である場合には他の前記演算装置を前記コピー先として決定する
ことを特徴とする情報処理システム。
【請求項9】
行列積を用いて最適化問題の最適解探索を行う情報処理システムが行う情報処理方法であって、前記最適化問題を、第１の状態変数ベクトルと所定行列と第２の状態変数ベクトルとの積を含んだ目的関数で表される二次計画問題に変換する変換処理と、
前記所定行列と前記第１の状態変数ベクトルとを乗算して第１ベクトルを計算し、該第１ベクトルに基づくパラメータを含んだ確率分布に基づいて前記第２の状態変数ベクトルの確率的更新を行い、前記所定行列と前記第２の状態変数ベクトルとを乗算して第２ベクトルを計算し、該第２ベクトルに基づくパラメータを含んだ確率分布に基づいて前記第１の状態変数ベクトルの確率的更新を行う処理を繰り返すことで、前記第１の状態変数ベクトル及び前記第２の状態変数ベクトルの確率的更新を行う状態更新処理と、
前記第１の状態変数ベクトルと前記第２ベクトルとの内積によって前記積を計算することで、前記目的関数の関数値を算出する関数値算出処理と
を含んだことを特徴とする情報処理方法。
【請求項10】
請求項９に記載の情報処理方法であって、前記状態更新処理及び前記関数値算出処理は、
前記第１の状態変数ベクトル及び前記第２の状態変数ベクトルに対応する複数のレプリカ毎に独立かつ並行に実行される
ことを特徴とする情報処理方法。
【請求項11】
請求項１０に記載の情報処理方法であって、前記情報処理システムは、複数の演算装置を有し、
前記レプリカ毎の前記状態更新処理及び前記関数値算出処理は、前記複数の演算装置によって実行される
ことを特徴とする情報処理方法。
【請求項12】
請求項１１に記載の情報処理方法であって、前記関数値算出処理によって前記レプリカ毎の前記関数値を算出し、
前記レプリカ毎の前記関数値及び逆温度の差分に基づいて前記レプリカ毎の重みを算出し、該重みに基づいて、前記複数のレプリカから選択するレプリカの選択数を決定するリサンプリング処理
を含んだことを特徴とする情報処理方法。
【請求項13】
請求項１２に記載の情報処理方法であって、前記リサンプリング処理において、前記選択数に応じて、前記選択するレプリカが配置される前記演算装置に該レプリカのコピーを配置可能である場合には該演算装置が該レプリカのコピーを配置するコピー先として決定され、配置不可能である場合には他の前記演算装置が前記コピー先として決定される
ことを特徴とする情報処理方法。
【請求項14】
請求項１～８の何れか１項に記載の情報処理システムとしてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、情報処理システム、情報処理方法、及び情報処理プログラムに関する。
【背景技術】
【0002】
特許文献１には、「イジングモデルの１つのスピンを３以上の状態で表現する値を記憶する第１のメモリセルと、１つのスピンに相互作用を及ぼす他のスピンからの相互作用を示す相互作用係数を記憶する第２のメモリセルと、他のスピンの状態を表現する値と前記相互作用係数を定数または変数として持つ関数に基づいて、１つのスピンの次状態を決定する論理回路と、を有する単位ユニットを複数備える半導体装置」が開示されている。
【0003】
また、特許文献２には、任意の結合を持つイジングモデルに対して、マルコフ連鎖モンテカルロ法の要求する理論的背景を満たしつつ、全スピンを同時に確率的更新して最適解探索を実現する方法が開示されている。
【先行技術文献】
【特許文献】
【0004】
【文献】特開２０１６－５１３１４号公報
【文献】国際公開第２０１９／２１６２７７号公報
【非特許文献】
【0005】
【文献】Okuyama, T., Sonobe, T., Kawarabayashi, K. I., &Yamaoka, M. (2019). Binary optimization by momentum annealing. Physical ReviewE, 100(1), 012111.
【文献】Botev, Z. I. (2017). The normal law under linearrestrictions: simulation and estimation via minimax tilting. Journal of theRoyal Statistical Society: Series B (Statistical Methodology), 79(1), 125-148.
【文献】Neal, R. M. (1998). Suppressing random walks in Markovchain Monte Carlo using ordered overrelaxation. In Learning in graphical models(pp. 205-228). Springer, Dordrecht.
【発明の概要】
【発明が解決しようとする課題】
【0006】
物理現象や社会現象の多くは相互作用モデルで表現可能である。相互作用モデルは、モデルを構成する複数のノードと、ノード間の相互作用、さらに必要であればノード毎に作用する係数で定義される。物理学や社会科学の分野においては、イジングモデルを始めとする種々のモデルが提案されているが、いずれも相互作用モデルの一形態として解釈することができる。
【0007】
この相互作用モデルに関係づけられた指標を最小化又は最大化するノード状態を求めることが社会課題の解決において重要である。例えば、ソーシャルネットワークのクリークを検知する問題や、金融分野のポートフォリオ最適化問題が挙げられる。これらは、オペレーションズ・リサーチの分野では、無制約二値二次計画問題や混合二値二次計画問題に大別される。
【0008】
本発明は上述の背景に鑑みてなされたもので、相互作用モデルの基底状態探索を含む最適化問題の最適解探索を高速化する技術の提供を目的とする。
【課題を解決するための手段】
【0009】
上述した課題を解決するため、本発明の一態様では、行列積を用いて最適化問題の最適解探索を行う情報処理システムであって、記憶部と協働して処理を実行する処理部を有し、前記処理部は、前記最適化問題を、第１の状態変数ベクトルと所定行列と第２の状態変数ベクトルとの積を含んだ目的関数で表される二次計画問題に変換する変換処理と、前記所定行列と前記第１の状態変数ベクトルとを乗算して第１ベクトルを計算し、該第１ベクトルに基づくパラメータを含んだ確率分布に基づいて前記第２の状態変数ベクトルの確率的更新を行い、前記所定行列と前記第２の状態変数ベクトルとを乗算して第２ベクトルを計算し、該第２ベクトルに基づくパラメータを含んだ確率分布に基づいて前記第１の状態変数ベクトルの確率的更新を行う処理を繰り返すことで、前記第１の状態変数ベクトル及び前記第２の状態変数ベクトルの確率的更新を行う状態更新処理と、前記第１の状態変数ベクトルと前記第２ベクトルとの内積によって前記積を計算することで、前記目的関数の関数値を算出する関数値算出処理とを実行することを特徴とする。
【発明の効果】
【0010】
本発明によれば、相互作用モデルの基底状態探索を含む最適化問題の最適解探索を高速化できる。
上記した以外の課題、構成及び効果は、以下の発明を実施するための形態の説明により明らかにされる。
上記した以外の課題、構成及び効果は、以下の発明を実施するための形態の説明により明らかにされる。
【図面の簡単な説明】
【0011】
【図1】最適化問題の変数配列及び目的関数値の関係を例示する概念図である。
【図2】最適化問題の目的関数の各変数間の関係を例示する図である。
【図3】最適化問題の目的関数の各変数間の関係を例示する図である。
【図4】複数のレプリカにおける最適解探索処理を例示する図である。
【図5】実施形態に係る情報処理装置を実現するコンピュータのハードウェアを例示する図である。
【図6】実施形態に係る演算装置を構成する演算回路を例示する機能ブロック図である。
【図7】実施形態に係る情報処理装置の構成を例示する機能ブロック図である。
【図8】実施形態に係る情報処理装置が実行する最適解探索処理を例示するフローチャートである。
【図9】相互作用演算処理の詳細を例示するフローチャートである。
【図10】リサンプリングによるレプリカ数の推移のシミュレーション結果を例示す図である。
【図11】実施形態による計算速度を比較例と比較したシミュレーション結果を例示する図である。
【発明を実施するための形態】
【0012】
以下、図面を参照して本発明の実施形態を説明する。実施形態は、本発明を説明するための例示であって、説明の明確化のため、適宜、省略及び簡略化がなされている。本発明は、他の種々の形態でも実施することが可能である。特に限定しない限り、各構成要素は単数でも複数でもよい。
【0013】
同一あるいは同様の機能を有する構成要素が複数ある場合には、同一の符号に異なる添字を付して説明する場合がある。また、これらの複数の構成要素を区別する必要がない場合には、添字を省略して説明する場合がある。
【0014】
実施形態において、プログラムを実行して行う処理について説明する場合がある。ここで、計算機は、プロセッサ（例えばＣＰＵ、ＧＰＵ）によりプログラムを実行し、記憶資源（例えばメモリ）やインターフェースデバイス（例えば通信ポート）等を用いながら、プログラムで定められた処理を行う。そのため、プログラムを実行して行う処理の主体を、プロセッサとしてもよい。同様に、プログラムを実行して行う処理の主体が、プロセッサを有するコントローラ、装置、システム、計算機、ノードであってもよい。プログラムを実行して行う処理の主体は、演算部であればよく、特定の処理を行う専用回路を含んでいてもよい。ここで、専用回路とは、例えばＦＰＧＡ（Field Programmable Gate Array）やＡＳＩＣ（Application Specific Integrated Circuit）、ＣＰＬＤ（Complex Programmable Logic Device）等である。
【0015】
プログラムは、プログラムソースから計算機にインストールされてもよい。プログラムソースは、例えば、プログラム配布サーバ又は計算機が読み取り可能な記憶メディアであってもよい。プログラムソースがプログラム配布サーバの場合、プログラム配布サーバはプロセッサと配布対象のプログラムを記憶する記憶資源を含み、プログラム配布サーバのプロセッサが配布対象のプログラムを他の計算機に配布してもよい。また、実施形態において、２以上のプログラムが１つのプログラムとして実現されてもよいし、１つのプログラムが２以上のプログラムとして実現されてもよい。
【0016】
（実施形態の理論的背景）
相互作用モデルに関連付けられた目的関数を最適化する最適化問題の変数をｓ１，…，ｓＮのＮ個とし、各変数ｓｉの定義域Ｄｉを２値｛－１，＋１｝又は連続値［－１，＋１]のいずれかであるとする。各変数ｓｉの定義域Ｄｉがどちらであるかは問題毎に決定される。そして、最適化問題の目的関数Ｈは式１で表される。すなわち、目的関数Ｈが変数ｓの２次式で表される。
相互作用モデルに関連付けられた目的関数を最適化する最適化問題の変数をｓ１，…，ｓＮのＮ個とし、各変数ｓｉの定義域Ｄｉを２値｛－１，＋１｝又は連続値［－１，＋１]のいずれかであるとする。各変数ｓｉの定義域Ｄｉがどちらであるかは問題毎に決定される。そして、最適化問題の目的関数Ｈは式１で表される。すなわち、目的関数Ｈが変数ｓの２次式で表される。
【0017】
【数1】
【0018】
【数2】
【0019】
【数3】
【0020】
【数4】
【0021】
もしすべてのｉに対してＤｉ＝｛－１，＋１｝ならば、この最適化問題はイジングモデルの基底状態探索問題と呼ばれる組合せ最適化問題である。本実施形態では、イジングモデルの基底状態探索を含む最適化問題において、マルコフ連鎖モンテカルロ法（以降、ＭＣＭＣ（Markov Chain Monte Carlo methods）と称する）を活用したアルゴリズムで最適解もしくは近似解を探索する。
【0022】
図１は、最適化問題の変数配列及び目的関数値の関係を例示する概念図であり、変数配列に対する目的関数値のランドスケープを表す。グラフの横軸は変数配列ｓ、縦軸は目的関数Ｈ（ｓ）の値である。ＭＣＭＣは現在の状態ｓから、状態ｓの近傍のある状態ｓ’への確率的な遷移を繰り返す。状態ｓから状態ｓ’に遷移する確率を、遷移確率Ｐ（ｓ，ｓ’）と称する。遷移確率Ｐの例としてメトロポリス法（Metropolis method）や熱浴法（heat-bath algorithm）が挙げられる。
【0023】
遷移確率は温度と呼ばれるパラメータを有し、これは状態間の遷移のし易さを表す。温度を大きな値から徐々に減少させつつＭＣＭＣを実行するとき、目的関数値が最も低い状態（最低エネルギーの状態、図１のＸ）に漸近的に収束する。これを利用して最小化問題の最適解又は近似解を求める手法が、シミュレーティッド・アニーリング（以下、ＳＡ（Simulated Annealing）と称する）や非特許文献１で提案されたモメンタム・アニーリング（以下、ＭＡ（Momentum Annealing）と称する）である。
【0024】
式４に示す最小化問題を解くにあたり、代わりに式５の最小化問題を解くことを考える。ただし集合Ｓｒｅｌａｘｅｄ＝｛ｓ|ｓi∈［－１，＋１]｝である。
【0025】
【数5】
【0026】
式５の最適解をｓ＊＝［ｓ１
＊，…，ｓＮ
＊］と表す。証明は割愛するが、式６で求まるｓ＋＝［ｓ１
＋，…，ｓＮ
＋］は式４の最適解の一つとなる。本実施形態の目標は式２の最適解探索であるが、式５の最適解ｓ＊を求解後に式６の変換を得ても、所望の解ｓ＋を得られるということである。ただし、関数ｓｇｎは引数が０以上ならば＋１、それ以外ならば－１を返す関数である。
【0027】
【数6】
【0028】
ここで、Ｎ次元ベクトルｖ＝［ｖ１，…，ｖＮ］を導入して、式７に示す関数Ｈ’を定義する。
【0029】
【数7】
【0030】
ただし、関数Ｖ（ｖ）は式８に記す定義の通りである。
【0031】
【数8】
【0032】
行列Ｗ＝ｄｉａｇ（ｗ１，…，ｗＮ）は任意の対角行列で、ｖiは［－１，＋１］を動く実数である。式５に示す最小化問題の代わりにＨ’(ｓ、ｖ)の最小化問題である式９を導入する。
【0033】
【数9】
【0034】
２つのＮ次元ベクトルｘ＝s＋ｖ、ｙ＝s－ｖを定義する。本来解きたい最適化問題の目的関数はＨのみだが、ここにＶという関数を導入することで、ＭＣＭＣで並列更新可能な関数を新たに得られるようにしている。すると、関数Ｈ’は式１０と書き直せる。
【0035】
【数10】
【0036】
つまり、式５の最小化問題は式１１の最小化問題と言い換えられる。
【0037】
【数11】
【0038】
式１１の最適解をｘ＊、ｙ＊と表すと、s＊＝（ｘ＊＋ｙ＊）／２なる等式が成り立つ。これらの議論はＷが零行列でも成り立つ。
【0039】
以上より、式２で表す最適化問題の最適解は、式１１に示す制約付き二次計画問題の解から求められる。この解を求めるために、ＭＣＭＣを活用する。
【0040】
図２は、最適化問題の目的関数の各変数間の関係を例示する図であり、式１１における関数Ｇの各変数どうしの関係を示したグラフィカルモデルを表す。図２では、Ｎ＝６の場合を例示しているが、一般的なＮについても同様である。関数Ｇの各変数どうしの関係は、完全２部グラフで表すことができる。関数Ｇ内で変数ｘｉに乗ぜられる変数は、ｙ１，…，ｙＮとｘｉのみである。ＭＣＭＣは変数値を確率的に更新するとき、その変数に係わる変数の値を用いる。つまり、変数ｘ１の値を更新するときｙ１，…，ｙＮ及びｘ１を求め、それ以外の変数（ここではｘ２，…，ｘＮ）を参照しない。これは他の変数、例えばｘ２の値の更新でも同様である。ゆえに、変数配列ｙの値が一定ならば、配列ｘのそれぞれの値を独立に同時に確率的更新してもＭＣＭＣの理論的要請は破らない。
【0041】
同様に変数ｙｉに乗ぜられる変数も、ｘ１，…，ｘＮとｙｉのみである。ゆえに、変数配列ｘの値が一定の下で、配列ｙのそれぞれの値を独立に同時に確率的更新できる。
【0042】
以上より、「ｘ１，…，ｘＮの同時更新」と「ｙ１，…，ｙＮの同時更新」を繰り返す手続きで構成されたＭＣＭＣを実行することで、並列化による高速化という利点を享受しながら関数Ｇを最小化する配列ｘ、ｙを探索できる。
【0043】
本実施形態の議論では、行列Ｊに制約を設けていないことに注意されたい。たとえば行列Ｊの全要素が非零である場合にも、上記の議論が成り立つため、並列更新が可能である。
【0044】
図３は、最適化問題の目的関数の各変数間の関係を例示する図であり、全結合グラフの例である。図３では、Ｎ＝６の場合を例示しているが、一般的なＮについても同様である。一方で、原問題である式２の最小化問題に対して直接ＭＣＭＣを適用する場合、変数配列ｓの係わり方が図３に示すように全結合グラフで表現されるため、一度に一変数しか確率更新できず、逐次更新に限定される。
【0045】
ここからは、各変数に対する確率的更新の手続きを述べる。更新対象の変数をｘｉとする。変数ｙ１，…，ｙＮの値が一定下では、温度Ｔのボルツマン分布における変数ｘｉの存在確率ｐ（ｘｉ）は式１２を満たす。
【0046】
【数12】
【0047】
ただし、変数Ａiは式１３で求める値である。
【0048】
【数13】
【0049】
変数ｘｉとｙｉは｜ｘｉ｜＋｜ｙｉ｜≦２であるため、ｘｉを動かせる範囲は－（２－｜ｙｉ｜）以上（２－｜ｙｉ｜）以下である。よって、変数ｘｉは平均Ａi／ｗi、分散Ｔ／ｗiの正規分布で－（２－｜ｙｉ｜）以上（２－｜ｙｉ｜）以下を定義域とする切断正規分布を基に、ｘｉの次状態をサンプリングすればよい。この方法ではｘｉの現在の状態には依らずに次状態を決めるということである。ｙｉについても同様である。本明細書では、ｘとｙの変数を区別しない場合ｓと表記することがある。
【0050】
標準正規分布に従う乱数はBox-Muller法で生成可能である。ここでは定義域が限定されるため、非特許文献２で示されたアルゴリズムを用いればよい。
【0051】
最適解探索は、温度０における平衡状態からのサンプリングと見なせる。そのため、良質な解探索の実現には、平衡状態への短時間での収束が好ましい。平衡状態への収束性を高めるため、ＭＣＭＣでは様々な技術が提案されており、これらの活用も可能である。たとえば、非特許文献３は過剰緩和法を提案している。これは次状態の候補として、温度Ｔのボルツマン分布から１つだけではなく、Ｋ個の状態をサンプリングする。そして、サンプリングしたＫ個の状態に加えて、現在の状態の計（Ｋ＋１）個の状態を並び替えてｘｃ
０≦…≦ｘｃ
r＝ｘｉ≦ｘｃ
Ｋと表す。つまり、現在の状態は（Ｋ＋１）個の値のうち、小さい方から（ｒ＋１）番目ということである。そしてｘｃ
Ｋ＋１－ｒを次状態に採用する。この方法では、次状態がｘｉの現在の状態に依存する。
【0052】
さて、上述のように、最適解探索は、相互作用モデルの状態を表す変数ｘ、ｙの確率的遷移を繰り返す確率的探索プロセスによって行われる。しかし、確率的遷移の結果、局所解に陥り、多くの試行回数を経ても最適解もしくは近似解が求まらないケースがある。そこで、本実施形態では、同一モデルの異なる状態にそれぞれ対応する複数のレプリカを作成し、レプリカ毎に確率的探索プロセスを並列かつ独立に実行する。そして、目的関数値に基づいて状態の評価を行い、評価が高いレプリカを優先的に選択して複製するリサンプリングを行うことで、局所解に陥ることなく大局的な状態更新を行って、最適解探索の高速化を図る。
【0053】
【0054】
例えば温度Ｔ＝Ｔ０（時間ステップ０）でそれぞれの初期状態がセットされた複数のノード、例えば３個のノード＃１、＃２、＃３に配置されるレプリカ１－１、１－２、１－３を用いて、独立かつ並行して最適解探索が行われる。温度Ｔ＝Ｔｔ（時間ステップｔ）になると、それまで独立かつ並行して最適解探索が行われたレプリカ１－１、１－２、１－３の目的関数値（式１１の関数Ｇ）がそれぞれ算出される。Ｒ個のレプリカのうちのｊ番目のレプリカ１－ｊ（ｊ＝１，・・・，Ｒ）の目的関数を関数Ｇｊとする。
【0055】
続いて、レプリカ１－ｊ（ｊ＝１，・・・，Ｒ）の重みαｊを、式１４のように算出する。
【0056】
【数14】
【0057】
ただし、Δβｔは、リサンプリングを行う際の時間ステップｔにおける温度Ｔｔの逆温度βｔ（＝１／Ｔｔ）の差分（βｔ－βｔ－１）である。このレプリカ１－ｊの重みαｊから、リサンプリング後の各レプリカ１－ｊの出現頻度を表す出現確率ｑｊを式１５のように算出する。
【0058】
【数15】
【0059】
ただし、式１５の右辺の分母は、リサンプリング前の全てのレプリカ１－ｊの重みαｊの和を表す。すなわち、式１５は、重みαｊを正規化する。
【0060】
そして、リサンプリング後のレプリカとして、各レプリカ１－ｊの選択数Ｎｊが式１６のように決定される。選択数Ｎｊとは、次の時間ステップの最適解探索で用いるレプリカ１－ｊの個数である。Ｎｊの小数点以下は適宜端数処理される。Ｒはレプリカの総数を表し、リサンプリングの前後で一定とするが、リサンプリング後に増加又は減少するとしてもよい。
【0061】
【数16】
【0062】
そして、選択数Ｎｊに応じてリサンプリング後のレプリカのコピー先が決定され、必要に応じて各レプリカ１－ｊの情報がコピーされる。
【0063】
図４の例示では、Ｒ＝３である。レプリカ１－１はＲ×ｑ１＝０．６、レプリカ１－２はＲ×ｑ２＝２．１、レプリカ１－３はＲ×ｑ３＝０．０３であるので、四捨五入によってレプリカ１－１が１個、レプリカ１－２が２個、レプリカ１－３が０個とそれぞれ選択される。すなわち、リサンプリング後のレプリカは、１個のレプリカ１－１及び２個のレプリカ１－２である。レプリカ１－１及び１個目のレプリカ１－２は、リサンプリング前から引き続いて用いられる。そして、２個目のレプリカ１－２は、ノード＃２には配置できないため、レプリカ１－３が選択されないノード＃３にレプリカ１－２の情報が転送される。リサンプリング後、１個のレプリカ１－１及び２個のレプリカ１－２のそれぞれについて確率的探索プロセスが並列かつ独立に実行される。
【0064】
最適解探索においてレプリカが並列に実行され、所定温度まで低下する等の実行終了条件が充足されると、複数のレプリカのうちの関数Ｇの値が最小であるレプリカの状態を表す変数に基づいて最適解もしくは近似解が出力される。
【0065】
【0066】
（情報処理装置１０のハードウェア構成）
図５は、実施形態に係る情報処理装置１０を実現するコンピュータのハードウェアを例示する図である。情報処理装置１０は、ハードウェアとして、プロセッサ１１、主記憶装置１２、補助記憶装置１３、入力装置１４、出力装置１５、通信装置１６、１又は複数の演算装置１７、及びこれらの装置を通信可能に接続するシステムバス１８を有する。情報処理装置１０は、例えば、その一部又は全部がクラウドシステム（Cloud System）により提供されるクラウドサーバ（Cloud Server）のような仮想的な情報処理資源を用いて実現されるものであってもよい。また情報処理装置１０は、例えば、互いに協調して動作する、通信可能に接続された複数の情報処理装置によって実現されるものであってもよい。
図５は、実施形態に係る情報処理装置１０を実現するコンピュータのハードウェアを例示する図である。情報処理装置１０は、ハードウェアとして、プロセッサ１１、主記憶装置１２、補助記憶装置１３、入力装置１４、出力装置１５、通信装置１６、１又は複数の演算装置１７、及びこれらの装置を通信可能に接続するシステムバス１８を有する。情報処理装置１０は、例えば、その一部又は全部がクラウドシステム（Cloud System）により提供されるクラウドサーバ（Cloud Server）のような仮想的な情報処理資源を用いて実現されるものであってもよい。また情報処理装置１０は、例えば、互いに協調して動作する、通信可能に接続された複数の情報処理装置によって実現されるものであってもよい。
【0067】
プロセッサ１１は、例えば、ＣＰＵ（Central Processing Unit）やＭＰＵ（Micro Processing Unit）を用いて構成されている。主記憶装置１２は、プログラムやデータを記憶する装置である。例えば、ＲＯＭ（Read Only Memory）、ＳＲＡＭ（Static Random Access Memory）、ＮＶＲＡＭ（Non Volatile RAM）、マスクＲＯＭ（Mask Read Only Memory）、ＰＲＯＭ（Programmable ROM）等）、ＲＡＭ（Random Access Memory）、ＤＲＡＭ（Dynamic Random Access Memory）等である。補助記憶装置１３は、ＨＤＤ（Hard Disk Drive）、フラッシュメモリ（Flash Memory）、ＳＳＤ（Solid State Drive）、光学式記憶装置（ＣＤ（Compact Disc）、ＤＶＤ(Digital Versatile Disc)等）等である。補助記憶装置１３に格納されているプログラムやデータは、随時、主記憶装置１２に読み込まれる。
【0068】
入力装置１４は、ユーザから情報の入力を受け付けるユーザインタフェースであり、例えば、キーボード、マウス、カードリーダ、タッチパネル等である。出力装置１５は、ユーザに情報を提供するユーザインタフェースであり、例えば、各種情報を可視化する表示装置（ＬＣＤ（Liquid Crystal Display）、グラフィックカード等）や音声出力装置（スピーカ）、印字装置等である。通信装置１６は、他の装置と通信する通信インタフェースであり、例えば、ＮＩＣ（Network Interface Card）、無線通信モジュール、ＵＳＢ（Universal Serial Interface）モジュール、シリアル通信モジュール等である。
【0069】
演算装置１７は、組合せ最適化問題の最適解探索に関する処理を実行する装置である。演算装置１７は、例えば、ＧＰＵ（Graphics Processing Unit）のように、情報処理装置１０に装着する拡張カードの形態を取るものであってもよい。演算装置１７は、例えば、ＣＭＯＳ（Complementary Metal Oxide Semiconductor）回路、ＦＰＧＡ）、ＡＳＩＣ等のハードウェアによって構成される。
【0070】
演算装置１７は、制御装置、記憶装置、システムバス１８に接続するためのインタフェース等を含み、システムバス１８を介してプロセッサ１１との間でコマンドや情報の送受を行う。演算装置１７は、例えば、通信線を介して他の演算装置１７と通信可能に接続され、他の演算装置１７と協調して動作するものであってもよい。演算装置１７により実現される機能を、例えば、プロセッサ（ＣＰＵ、ＧＰＵ等）にプログラムを実行させることにより実現してもよい。
【0071】
情報処理装置１０において、プログラムが補助記憶装置１３から読み出されて、プロセッサ１１及び主記憶装置１２の協働により実行されることにより、情報処理装置１０の機能が実現される。あるいは、情報処理装置１０の機能を実現するためのプログラムは、通信装置１６を介した通信により外部のコンピュータから取得されてもよい。あるいは、情報処理装置１０の機能を実現するためのプログラムは、可搬型の記録媒体（光学ディスク、磁気ディスク、光磁気ディスク、半導体記憶媒体等）に記録され、媒体読取装置により読み出されてもよい。
【0072】
また、情報処理装置１０は、複数の装置が協働して処理を実行する情報処理システムであってもよい。情報処理システムを実現するプログラムは、各プログラムによって各装置に各機能を実現させることで、情報処理装置１０と同様の機能を実現する。
【0073】
（演算回路５００の構成）
図６は、実施形態に係る演算装置１７を構成する演算回路５００を例示する機能ブロック図であり、演算装置１７の動作原理を示す。各演算装置１７は、１又は複数の演算回路５００を含んで構成される。各演算回路５００は、１つのレプリカに対応する。本実施形態では、演算回路５００は、例えば後述の差分計算ブロック５１４、サンプリングブロック５１５、次状態決定ブロック５１６等が、ＧＰＵやベクトル型計算機などに搭載された多数の演算器等を用いるソフトウェアで実装されるとする。しかし、これに限らず、ハードウェアで実装されてもよい。演算回路５００は、変数配列ｘ１，…，ｘＮ又は変数配列ｙ１，…，ｙＮを温度Ｔにおけるボルツマン分布（式１２）からサンプリングする機能を実現する。
図６は、実施形態に係る演算装置１７を構成する演算回路５００を例示する機能ブロック図であり、演算装置１７の動作原理を示す。各演算装置１７は、１又は複数の演算回路５００を含んで構成される。各演算回路５００は、１つのレプリカに対応する。本実施形態では、演算回路５００は、例えば後述の差分計算ブロック５１４、サンプリングブロック５１５、次状態決定ブロック５１６等が、ＧＰＵやベクトル型計算機などに搭載された多数の演算器等を用いるソフトウェアで実装されるとする。しかし、これに限らず、ハードウェアで実装されてもよい。演算回路５００は、変数配列ｘ１，…，ｘＮ又は変数配列ｙ１，…，ｙＮを温度Ｔにおけるボルツマン分布（式１２）からサンプリングする機能を実現する。
【0074】
演算回路５００は、変数メモリ５１１、非線形係数メモリ５１２、線形係数メモリ５１３、差分計算ブロック５１４、サンプリングブロック５１５、及び次状態決定ブロック５１６を含む。
【0075】
各演算回路５００の変数メモリ５１１には、前述した変数ｘ１，…，ｘＮ及びｙ１，…，ｙＮを示す情報が格納される（図２参照）。
【0076】
非線形係数メモリ５１２には、行列Ｊを表す情報が格納される。行列Ｊは一般に対称行列であり、この対称性を用いて非線形係数メモリ５１２の使用量を削減することができる。線形係数メモリ５１３には、ベクトルｈを表す情報が格納される。
【0077】
演算装置１７には、制御信号ＥＮ、重み信号ＳＷ、及び温度信号ＴＥが入力される。
【0078】
制御信号ＥＮは、Ｈ（high）とＬ（low）の値を周期的に繰り返す信号で、変数配列ｘとｙいずれを更新しているかを表す。たとえば、ＥＮがＨのときは変数配列ｘを更新、Ｌのときはｙを更新と定める。この制御信号ＥＮにより、変数ｘ１，…，ｘＮを同時に更新し、また変数ｙ１，…，ｙＮを同時に更新する。図５では簡略化のため制御信号ＥＮはサンプリングブロック５１５のみに入力しているが、変数メモリ５１１等の本信号を必要とする他の箇所に対しても同様に入力される。
【0079】
重み信号ＳＷは、対角行列Ｗの対角成分を表すＮ要素のベクトルを表す信号である。
【0080】
差分計算ブロック５１４には、非線形係数メモリ５１２に格納されている行列Ｊの値、線形係数メモリ５１３に格納されているベクトルｈ、重み信号ＳＷ、及び変数メモリ５１１に格納されている変数ｓ（ｘ又はｙ）が入力される。差分計算ブロック５１４は、制御信号ＥＮがＨのとき（Ｊ＋ｄｉａｇ（ｗ１，・・・，ｗＮ））ｙ＋ｈ、ＥＮがＬのとき（Ｊ＋ｄｉａｇ（ｗ１，・・・，ｗＮ））ｘ＋ｈを出力する。この出力値は前述のＡｉに相当する。
【0081】
サンプリングブロック５１５は、差分計算ブロック５１４の出力と重み信号ＳＷ、温度パラメータの値を保持する温度信号ＴＥ、制御信号ＥＮ、及び他の変数の値を受けとる。そしてｉ番目の要素として、制御信号ＥＮがＨのとき－（２－｜ｙｉ｜）以上（２－｜ｙｉ｜）以下、ＥＮがＬのとき－（２－｜ｘｉ｜）以上（２－｜ｘｉ｜）以下を定義域とする、式１２で表される切断正規分布からランダムにサンプリングして出力する。
【0082】
次状態決定ブロック５１６は、サンプリングブロック５１５から出力される１又は複数の値を基に、変数の次状態を決定する。もし、ＭＣＭＣの更新則として単なる熱浴法に定めたならば、次状態決定ブロック５１６はサンプリングブロック５１５の出力値を１つだけ受け取り、それをそのまま変数メモリ５１１に書き込めばよい。また、ＭＣＭＣの更新則として公知の過剰緩和法を用いるならば、次状態決定ブロック５１６はサンプリングブロック５１５から複数の値、そして変数メモリ５１１から更新対象の変数の現在値を受け取り、過剰緩和法に従って１つ選択して、変数メモリ５１１に書き込む。周知のように、過剰緩和法では、直前の状態との相関が負となるように次の状態を決める。
【0083】
（情報処理装置１０の機能構成）
図７は、実施形態に係る情報処理装置１０の構成を例示する機能ブロック図であり、情報処理装置１０が備える主な機能（ソフトウェア構成）を示す。情報処理装置１０は、処理部１０１、記憶部１０２、及び、１又は複数の演算装置１７を有する。
図７は、実施形態に係る情報処理装置１０の構成を例示する機能ブロック図であり、情報処理装置１０が備える主な機能（ソフトウェア構成）を示す。情報処理装置１０は、処理部１０１、記憶部１０２、及び、１又は複数の演算装置１７を有する。
【0084】
処理部１０１は、モデル変換部１０１ａ、モデル係数設定部１０１ｂ、重み設定部１０１ｃ、変数値初期化部１０１ｄ、温度設定部１０１ｅ、相互作用演算実行部１０１ｆ、及び変数値読出部１０１ｇを有する。これらの機能は、プロセッサ１１が主記憶装置１２と協働して補助記憶装置１３に格納されているプログラムを読み出して実行することにより、もしくは、演算装置１７が備えるハードウェアにより実現される。なお、情報処理装置１０は、上記の機能に加えて、例えば、オペレーティングシステム、ファイルシステム、デバイスドライバ、ＤＢＭＳ（Data-Base Management System）等の他の機能を備えていてもよい。
【0085】
記憶部１０２は、主記憶装置１２又は補助記憶装置１３により実現され、問題データ１０２ａ、二次計画形式問題データ１０２ｂ、定義域データ１０２ｃ、及び演算装置制御プログラム１０２ｄを記憶する。問題データ１０２ａは、例えば、最適化問題等を公知の所定の記述形式で記述したデータである。問題データ１０２ａは、例えば、ユーザがユーザインタフェース（入力装置、出力装置、通信装置等）を介して設定される。
【0086】
二次計画形式問題データ１０２ｂは、モデル変換部１０１ａが、問題データ１０２ａを式４が示す二次計画問題のフォーマットに合致する形式のデータに変換することにより生成されるデータである。この変換にあたり、与えられた各変数の定義域は、定義域データ１０２ｃに書き込まれる。定義域データ１０２ｃは、例えば各変数が２値を取るか実数値を取るかを示している。演算装置制御プログラム１０２ｄは、相互作用演算実行部１０１ｆが演算装置１７を制御する際に実行する、もしくは相互作用演算実行部１０１ｆが個々の演算装置１７にロードして演算装置１７に実行させるプログラムである。
【0087】
モデル変換部１０１ａは、問題データ１０２ａを二次計画問題のフォーマットである二次計画形式問題データ１０２ｂに変換する。このために、式１から式１１を導出する機能を、ソフトウェアあるいはハードウェアとしてモデル変換部１０１ａに実装しておけばよい。モデル変換部１０１ａの機能は必ずしも情報処理装置１０に実装されていなくてもよく、情報処理装置１０が、他の情報処理装置等で生成された二次計画形式問題データ１０２ｂを入力装置１４や通信装置１６を介して取り込むようにしてもよい。
【0088】
モデル係数設定部１０１ｂ、重み設定部１０１ｃ、変数値初期化部１０１ｄ、及び温度設定部１０１ｅは、演算装置１７にＲ個の演算回路５００を生成する。そして、Ｒ個の演算回路５００に対して以下の処理を実行することで、Ｒ個のレプリカ１－ｊ（ｊ＝１，２，・・・，Ｒ）を作成する。Ｒ個の演算回路５００は、１又は複数の演算装置１７に配置される。
【0089】
モデル係数設定部１０１ｂは、二次計画形式問題データ１０２ｂに基づく、式４の行列Ｊを非線形係数メモリ５１２に、ベクトルｈを線形係数メモリ５１３に設定する。
【0090】
変数値初期化部１０１ｄは、演算回路５００の変数メモリ５１１に格納されている各変数の値を初期化する。変数値初期化部１０１ｄは、例えば、各変数の値を－１以上＋１以下から一様に、ランダムサンプリングして決めればよい。この際に、変数に関する制約である｜ｘｉ｜＋｜ｙｉ｜≦２を満たすよう注意しなければならない。また、このときの各変数の値は連続値で扱われていることに留意されたい。
【0091】
温度設定部１０１ｅは、相互作用演算実行部１０１ｆが最適解探索を行う際に用いる温度Ｔを設定する。
【0092】
相互作用演算実行部１０１ｆは、温度設定部１０１ｅにより設定された各温度Ｔについて、式１１で表す関数Ｇを最小化する変数配列ｘ及びｙを探索する演算（以下、相互作用演算処理と称する。）をレプリカ１－ｊ毎の各演算回路５００に実行させる。相互作用演算処理に際し、相互作用演算実行部１０１ｆは、例えば、温度Ｔを高いほうから低い方に向けて変化させる。
【0093】
具体的には、各レプリカ１－ｊにおける変数配列の確率的な更新は、式１２、式１３のように、変数ｘの更新の際に変数ｙが用いられ、変数ｙの更新の際に変数ｘが用いられる。具体的には、式１２、式１３から分かるように、変数ｘの更新の際には、行列Ｊと変数ｙの乗算結果が用いられる。同様に、変数ｙの更新の際には、行列Ｊと変数ｘの乗算結果が用いられる。
【0094】
また、相互作用演算実行部１０１ｆは、リサンプリングを所定周期で実行する。すなわち、相互作用演算実行部１０１ｆは、時間ステップｔの温度Ｔにおいて、各演算回路５００（ノード）に作成されたレプリカ１－ｊ毎のエネルギー（目的関数値（関数Ｇｊの値））を計算する。
【0095】
ここで、関数Ｇｊの行列Ｊと変数ｙの乗算（式１１の右辺第１項）は、計算負荷が大きい。しかし、この行列Ｊと変数ｙの乗算は、変数ｘの前回の更新ステップの際に行われており、乗算結果Ｊｙがワークエリアに格納されている。そこで、関数Ｇｊの値の算出の際に、行列Ｊと変数ｙの乗算を行う代わりに、ワークメモリに格納されている乗算結果Ｊｙを用いる。乗算結果Ｊｙはベクトルであるので、関数Ｇｊの右辺第１項の計算は、ベクトルｘＴとの内積計算に置き換わる。このように、変数ｘ、ｙの更新過程の計算結果を目的関数値の算出に用いることで、目的関数値の計算量を抑制し、最適化探索の計算速度の高速化を図ることができる。
【0096】
そして、相互作用演算実行部１０１ｆは、時間ステップｔにおけるレプリカ１－ｊ毎のエネルギーをもとに、式１４、式１５、式１６から、リサンプリング後の各レプリカ１－ｊの選択数Ｎｊを決定する。
【0097】
そして、相互作用演算実行部１０１ｆは、各レプリカ１－ｊの選択数Ｎｊのうち、リサンプリング後の増加分の各レプリカ１－ｊのコピー先を決定する。各レプリカ１－ｊのコピー先は、コピー元のレプリカ１－ｊと同一の演算装置１７（自ノード）の上限内で自ノードが優先的に決定される。自ノードの上限を超過する場合には、レプリカのコピー先は、それぞれの上限内で自ノードとは異なる他の演算装置１７（他ノード）がランダムに決定される。上限とは、各ノードのリソース上限、又は、リサンプリング前に各ノードに配置されていたレプリカの数、である。そして、コピー先へコピー元のレプリカの情報が転送される。
【0098】
このように、自ノード内を優先的にコピー先とし、レプリカの情報のコピーの際に通信が発生する他ノードへのコピーを極力回避することで、情報転送の待ち時間の発生を抑制し、計算速度の低下を防止できる。
【0099】
変数値読出部１０１ｇは、相互作用演算実行部１０１ｆによる最適解探索が終了すると、関数Ｇの値が最小のレプリカに対応する演算回路５００の変数メモリ５１１に格納されている変数配列ｘ及びｙを読み出す。ここで読み出される値は、式１１の解である。上述の議論に従って、Ｎ次元ベクトルｓ＊＝（ｘ＋ｙ）／２を計算する。そして定義域データ１０２ｃを読み出し、式６で得られるベクトルｓ＋を最終的な解として出力装置１５や通信装置１６に出力する。つまり、定義域データ１０２ｃにてｉ番目の定義域が｛－１，＋１｝と判明すればｓｇｎ（ｓ＊ｉ）、ｉ番目の定義域が［－１，＋１］ならばｓｉ自体を出力するということである。このようにして、定義された値域に応じた解が求められる。
【0100】
（最適解探索処理のフロー）
図８は、実施形態に係る情報処理装置１０が実行する最適解探索処理を例示するフローチャートである。最適解探索処理は、例えば、入力装置１４を介してユーザからの指示等を受け付けることにより開始される。
図８は、実施形態に係る情報処理装置１０が実行する最適解探索処理を例示するフローチャートである。最適解探索処理は、例えば、入力装置１４を介してユーザからの指示等を受け付けることにより開始される。
【0101】
Ｓ１１では、モデル変換部１０１ａは、問題データ１０２ａを二次計画形式問題データ１０２ｂに変換する。二次計画形式問題データ１０２ｂでは、例えば式１で表現される目的関数Ｈにおける行列Ｊ、ベクトルｈが任意の形式で表現される。記憶部１０２が既に二次計画形式問題データ１０２ｂを記憶している場合には、Ｓ１１が省略される。Ｓ１１の処理と、Ｓ１２以降の処理とは、夫々を異なる装置で実行するようにしてもよい。またＳ１１の処理と、Ｓ１２以降の処理とを異なるタイミングで実行するようにしてもよい（例えば、Ｓ１１の処理を事前に行っておくことが考えられる）。
【0102】
Ｓ１２～Ｓ１５の処理は、演算装置１７のＲ個の演算回路５００にて作成されるＲ個のレプリカについて実行される。
【0103】
Ｓ１２では、モデル係数設定部１０１ｂは、非線形係数メモリ５１２及び線形係数メモリ５１３に行列Ｊとベクトルｈの値を設定する。メモリの値は、ユーザインタフェース（例えば、入力装置１４、出力装置１５、通信装置１６等により実現される）を介してユーザが設定又は編集することができる。
【0104】
Ｓ１３では、重み設定部１０１ｃは、重み信号ＳＷの値を決定する。前述の式８の説明通り、最適解を探索する上で重み信号ＳＷは任意の値を取ることが許される。そのため、信号値は常に０としてもよい。この場合は計算の負荷を軽減することができる。また、特許文献２の式３～式５に示すように行列Ｊの固有値から決定してもよい。あるいは、行列Ｊの行和から決定してもよい。重み信号ＳＷの値算出の計算は、演算装置１７内又はプロセッサ１１で実行してもよい。あるいはユーザが自分で設定してもよい（Ｓ１３）。
【0105】
Ｓ１４では、変数値初期化部１０１ｄは、変数メモリ５１１に格納されている各変数の値を初期化する。変数メモリ５１１に格納する値は連続値である。先に述べたように初期値はランダムでよい。以上で、式１１を表現するパラメータが設定されたことになる。
【0106】
Ｓ１５では、温度設定部１０１ｅは、最適解探索にて使用する温度パラメータの系列Ｔγ（γ＝１，２，３，…）を初期設定する。なお、上記の添字γは設定される温度Ｔの種類を表す。温度Ｔの設定方法については、例えば特許文献１の方法を採用可能である。
【0107】
Ｓ１６では、相互作用演算実行部１０１ｆは、相互作用演算処理を実行する。相互作用演算処理の詳細は、図９を参照して後述する。
【0108】
Ｓ１７では、相互作用演算実行部１０１ｆは、停止条件が成立したか否か（例えば、温度Ｔが予め設定された最低温度に達したか否か）を判定する。相互作用演算実行部１０１ｆは、停止条件が成立したと判定した場合（Ｓ１７：Ｙｅｓ）、Ｓ１８に処理を移す一方、停止条件が成立していないと判定した場合（Ｓ１７：Ｎｏ）、Ｓ１６に処理を戻す。
【0109】
Ｓ１８では、変数値読出部１０１ｇは、関数Ｇｊの値が最小のレプリカ１－ｊの演算回路５００から変数メモリ５１１に格納されている変数の値を読み出す。また、変数値読出部１０１ｇは、定義域データ１０２ｃに格納されている二次計画形式問題データ１０２ｂの各変数の定義域を読み出す。そして、変数値読出部１０１ｇは、式６に基づいた変換を通じたベクトルを算出して、式２もしくは式４の解として出力する。以上で最適解探索処理は終了する。
【0110】
（相互作用演算処理の詳細フロー）
図９は、図８のＳ１６の相互作用演算処理の詳細を例示するフローチャートである。Ｓ１６１では、相互作用演算実行部１０１ｆは、各演算回路５００の演算により、各レプリカ１－ｊの変数配列の確率的な同時更新を実行する。
図９は、図８のＳ１６の相互作用演算処理の詳細を例示するフローチャートである。Ｓ１６１では、相互作用演算実行部１０１ｆは、各演算回路５００の演算により、各レプリカ１－ｊの変数配列の確率的な同時更新を実行する。
【0111】
Ｓ１６２では、相互作用演算実行部１０１ｆは、レプリカ１－ｊ毎の関数Ｇｊ（式１１）の値を算出する。上述したように、レプリカ１－ｊ毎の関数Ｇｊの目標関数値の算出の際に、変数ｘ、ｙの更新過程の計算結果を用いる。なお、関数Ｇｊの値の算出は、各レプリカ１－ｊの演算回路５００が行ってもよい。
【0112】
Ｓ１６３では、相互作用演算実行部１０１ｆは、各レプリカ１－ｊの選択数Ｎｊ及び増加分のレプリカのコピー先を決定する。
【0113】
Ｓ１６４では、相互作用演算実行部１０１ｆは、Ｓ１６３で決定したコピー先のノード（自ノード又は他ノード）へコピーするレプリカの情報を送信する。
【0114】
（レプリカ数の推移のシミュレーション結果）
図１０は、リサンプリングによるレプリカ数の推移のシミュレーション結果を例示す図である。図１０では、縦軸は時刻を表し、横軸はレプリカ数を表す。図１０では、最も高い温度の時刻０でそれぞれ１個であったレプリカ１－１、１－２、１－３、１－４、１－５、１－６、１－７、１－８、１－９、１－１０がリサンプリングの繰り返しを経て増減する様子を表している。図１０では、時刻０におけるオリジナルのレプリカが同一であるレプリカは、全て同一の系統としている。そのため、時刻１以降で同一の塗りつぶし又はハッチングのパターンが施されているレプリカは、同じ変数を持っているわけではないことに注意されたい。例えば最も低い温度の時刻５０のレプリカ１－６及び１－９のように、一部の目的関数値が優良な系統のレプリカが多数を占めることで、目的関数値がより迅速に最低状態に収束する可能性が高まり、最適化問題をより高速処理できるようになる。
図１０は、リサンプリングによるレプリカ数の推移のシミュレーション結果を例示す図である。図１０では、縦軸は時刻を表し、横軸はレプリカ数を表す。図１０では、最も高い温度の時刻０でそれぞれ１個であったレプリカ１－１、１－２、１－３、１－４、１－５、１－６、１－７、１－８、１－９、１－１０がリサンプリングの繰り返しを経て増減する様子を表している。図１０では、時刻０におけるオリジナルのレプリカが同一であるレプリカは、全て同一の系統としている。そのため、時刻１以降で同一の塗りつぶし又はハッチングのパターンが施されているレプリカは、同じ変数を持っているわけではないことに注意されたい。例えば最も低い温度の時刻５０のレプリカ１－６及び１－９のように、一部の目的関数値が優良な系統のレプリカが多数を占めることで、目的関数値がより迅速に最低状態に収束する可能性が高まり、最適化問題をより高速処理できるようになる。
【0115】
（実施形態を含む開示技術の効果）
上述の実施形態を含む開示技術では、最適化問題を、第１の状態変数ベクトルと所定行列と第２の状態変数ベクトルとの積を含んだ目的関数で表される二次計画問題に変換し、目的関数の関数値を算出する際、第１の状態変数ベクトルと第２ベクトルとの内積によって目的関を計算する。ここで、第２ベクトルは、所定行列と第１の状態変数ベクトルとを乗算して第１ベクトルを計算し、第１ベクトルに基づくパラメータを含んだ確率分布に基づいて第２の状態変数ベクトルの確率的更新を行い、所定行列と第２の状態変数ベクトルとを乗算して第２ベクトルを計算し、第２ベクトルに基づくパラメータを含んだ確率分布に基づいて第１の状態変数ベクトルの確率的更新を行う処理を繰り返すことで第１の状態変数ベクトル及び第２の状態変数ベクトルの確率的更新を行った際の第２ベクトルが記憶部に記憶されたものである。よって、最適解探索の反復計算で現れた計算過程の結果を再利用することで、低コストで目的関数値を計算でき、最適解探索の処理を高速化できる。
上述の実施形態を含む開示技術では、最適化問題を、第１の状態変数ベクトルと所定行列と第２の状態変数ベクトルとの積を含んだ目的関数で表される二次計画問題に変換し、目的関数の関数値を算出する際、第１の状態変数ベクトルと第２ベクトルとの内積によって目的関を計算する。ここで、第２ベクトルは、所定行列と第１の状態変数ベクトルとを乗算して第１ベクトルを計算し、第１ベクトルに基づくパラメータを含んだ確率分布に基づいて第２の状態変数ベクトルの確率的更新を行い、所定行列と第２の状態変数ベクトルとを乗算して第２ベクトルを計算し、第２ベクトルに基づくパラメータを含んだ確率分布に基づいて第１の状態変数ベクトルの確率的更新を行う処理を繰り返すことで第１の状態変数ベクトル及び第２の状態変数ベクトルの確率的更新を行った際の第２ベクトルが記憶部に記憶されたものである。よって、最適解探索の反復計算で現れた計算過程の結果を再利用することで、低コストで目的関数値を計算でき、最適解探索の処理を高速化できる。
【0116】
また、開示技術では、状態更新及び関数値算出を、第１の状態変数ベクトル及び第２の状態変数ベクトルに対応する複数のレプリカ毎に独立かつ並行に実行するので、個々の局所解にとらわれずに、大局的に状態更新を行うことができる。
【0117】
また、開示技術では、レプリカ毎の状態更新及び関数値算出を、複数の演算装置によって実行するので、負荷を分散しながら並行処理を行って処理の高速化を図ることができる。
【0118】
また、開示技術では、関数値算出処理によってレプリカ毎の関数値を算出し、レプリカ毎の関数値及び逆温度の差分に基づいてレプリカ毎の重みを算出し、重みに基づいて、複数のレプリカから選択するレプリカの選択数を決定するリサンプリングを行う。よって、最適解探索の反復計算で現れた計算過程の結果を再利用することで、リサンプリングを繰り返す毎に関数値を算出する負荷が大きい計算を、Ｏ（log2ｎ）のオーダの低コストで目的関数値を計算でき、最適解探索の処理を高速化できる。
【0119】
また、開示技術では、リサンプリングにおけるレプリカの選択数は、リサンプリング前のレプリカの合計数と同一であることから、限られたリソースを目的関数値が高いレプリカの系列に多く配分することで、効率的に最適解探索を高速化できる。
【0120】
また、開示技術では、リサンプリングにおけるレプリカの選択数は、リサンプリング前の複数のレプリカの合計数よりも多いことから、潤沢なリソースを用いて多くの系のレプリカによって最適解探索を行うことで、最適解探索をさらに高速化できる。
【0121】
また、開示技術では、リサンプリングにおけるレプリカの選択数は、リサンプリング前の複数のレプリカの合計数よりも少ないことから、限られたリソースを目的関数値が高いレプリカの系列に配分することで、効率的に最適解探索を高速化できる。
【0122】
また、開示技術では、リサンプリングにおけるレプリカの選択数に応じて、選択するレプリカが配置される演算装置にこのレプリカのコピーを配置可能である場合にはこの演算装置をこのレプリカのコピーを優先的に配置するコピー先として決定する。一方、配置不可能である場合には他の演算装置をコピー先として決定する。よって、リサンプリングの際のレプリカの情報の転送（コピー）の際に、ノード間通信の発生を抑制するので、通信時間を抑制し、最適解探索を高速化できる。
【実施例】
【0123】
上述の実施形態の実施例として、実施例１：１ＧＰＵで実装（演算装置１７（ＧＰＵ）を１つ実装）、及び、実施例２：４ＧＰＵで実装（演算装置１７（ＧＰＵ）を４つ実装）、の場合の計算速度のシミュレーション結果を以下に示す。
【0124】
図１１は、実施形態による計算速度を比較例と比較したシミュレーション結果を例示する図である。図１１では、横軸は１回あたりのアニーリングの計算時間を表し、縦軸は目的関数値を表す。何れも乱数シードを変えて、合計１００回ずつの試行とし、各時間における目的関数値の平均値、第１四分位、及び第３四分をプロットした。破線が平均値を示し、破線に重なるグレー線の下端の外郭線が１四分位を示し、上端の外郭線が第３四分位を示す。
【0125】
表１に示すように、比較例１では、５０ノードでのＳＡ（Simulated Annealing）を２回行うことで合計１００回の試行とした。比較例２では、ＰＭＣ（Parallel Markov Chains）なしのＭＡ（Momentum Annealing）を１ＧＰＵで１００回試行した。比較例３では、ＰＭＣなしのＭＡを４ＧＰＵで１００回試行した。実施例１では、ＰＭＣありのＭＡを１ＧＰＵで１００回試行した。実施例２では、ＰＭＣありのＭＡを４ＧＰＵで１００回試行した。
【0126】
【表1】
【0127】
図１１から分かるように、実施例１及び２は、比較例１～３の何れと比較しても最低エネルギー状態（目的関数値＝－４．４０×１０５）付近により速く到達した。また、実施例１は、同じく１ＧＰＵ実装の比較例２と比較して、最低エネルギー状態付近により速く到達した。また、実施例２は、同じく４ＧＰＵ実装の比較例３と比較して、最低エネルギー状態付近により速く到達した。また、実施例２は、実施例１と比較して、最低エネルギー状態付近により速く到達した。
【0128】
以上のシミュレーション結果から、同じＭＡでも１ＧＰＵ実装よりは４ＧＰＵ実装の方が計算速度の面で有利であり、さらに本実施形態のＰＭＣありＭＡでは、ＧＰＵの実装数が多いほど計算速度面で有利であることが分かった。
【0129】
（その他の実施形態）
演算回路５００は、既に述べた最適化問題を解く計算を実行する機能を備える限り、ソフトウェアで構成してもよいし、ハードウェアで構成してもよい。具体的には、アニーリング方式において電子回路（デジタル回路など）で実装するハードウェアだけでなく、超伝導回路などで実装する方式でもよい。また、アニーリング方式以外にてイジングモデルを実現するハードウェアでもよい。例えばレーザーネットワーク方式（光パラメトリック発振）、量子ニューラルネットワークなどが知られている。また、一部の考え方が異なるものの、イジングモデルで行う計算をアダマールゲート、回転ゲート、制御ＮＯＴゲートといったゲートで置き換えた量子ゲート方式も、本実施形態の構成として採用することができる。
演算回路５００は、既に述べた最適化問題を解く計算を実行する機能を備える限り、ソフトウェアで構成してもよいし、ハードウェアで構成してもよい。具体的には、アニーリング方式において電子回路（デジタル回路など）で実装するハードウェアだけでなく、超伝導回路などで実装する方式でもよい。また、アニーリング方式以外にてイジングモデルを実現するハードウェアでもよい。例えばレーザーネットワーク方式（光パラメトリック発振）、量子ニューラルネットワークなどが知られている。また、一部の考え方が異なるものの、イジングモデルで行う計算をアダマールゲート、回転ゲート、制御ＮＯＴゲートといったゲートで置き換えた量子ゲート方式も、本実施形態の構成として採用することができる。
【0130】
演算回路５００の一つの実装例として、ＣＭＯＳ（Complementary Metal-Oxide Semiconductor）集積回路や、ＦＰＧＡ上の論理回路として実装することができる。例えば、特許文献１に開示されているように、ＳＲＡＭ（Static Random Access Memory）の技術を適用したユニットを多数配置し、各ユニットに変数を格納するメモリと変数を更新するための回路を配置してもよい。
【0131】
本発明は上述の実施形態に限定されるものではなく、様々な変形例を含む。例えば、上記した実施形態は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、矛盾しない限りにおいて、ある実施形態の構成の一部を他の実施形態の構成で置き換え、ある実施形態の構成に他の実施形態の構成を加えることも可能である。また、各実施形態の構成の一部について、構成の追加、削除、置換、統合、又は分散をすることが可能である。また、実施形態で示した構成及び処理は、処理効率又は実装効率に基づいて適宜分散、統合、又は入れ替えることが可能である。
【符号の説明】
【0132】
１０：情報処理装置、１１：プロセッサ、１２：主記憶装置、１３：補助記憶装置、１７：演算装置、５１１：変数メモリ、５１２：非線形係数メモリ、５１３：線形係数メモリ、５１４：差分計算ブロック、５１５：サンプリングブロック、５１６：次状態決定ブロック、１０１：処理部、１０１ａ：モデル変換部、１０１ｂ：モデル係数設定部、１０１ｃ：重み設定部、１０１ｄ：変数値初期化部、１０１ｅ：温度設定部、１０１ｆ：相互作用演算実行部、１０１ｇ：変数値読出部、１０２：記憶部、１０２ａ：問題データ、１０２ｂ：二次計画形式問題データ、１０２ｃ：定義域データ、１０２ｄ：演算装置制御プログラム
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
