知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024164601
(43)【公開日】2024-11-27
(54)【発明の名称】情報処理システム、情報処理方法及びプログラム
(51)【国際特許分類】
G06N 99/00 20190101AFI20241120BHJP
【ＦＩ】
G06N99/00 180
【審査請求】未請求
【請求項の数】7
【出願形態】ＯＬ
(21)【出願番号】P 2023080206
(22)【出願日】2023-05-15
(71)【出願人】
【識別番号】000005223
【氏名又は名称】富士通株式会社
(74)【代理人】
【識別番号】110002918
【氏名又は名称】弁理士法人扶桑国際特許事務所
(72)【発明者】
【氏名】穴澤 俊久
(72)【発明者】
【氏名】下川 聡
(57)【要約】
【課題】制約係数の決定にかかる時間を短縮する。
【解決手段】制約条件を含む組合せ最適化問題を、コスト項と制約条件に対応する制約項との和で表されるイジングモデルのエネルギー関数に基づいて計算する情報処理システム１０における処理部１２が、組合せ最適化問題の問題情報を取得し、問題情報に基づいて、制約条件を満たす解から制約条件を満たさない解への遷移が生じた場合のコスト項の変化量の期待値を計算し、期待値に基づいて制約項の係数を決定し、決定した係数を用いてエネルギー関数の情報を出力する。
【選択図】図１
【解決手段】制約条件を含む組合せ最適化問題を、コスト項と制約条件に対応する制約項との和で表されるイジングモデルのエネルギー関数に基づいて計算する情報処理システム１０における処理部１２が、組合せ最適化問題の問題情報を取得し、問題情報に基づいて、制約条件を満たす解から制約条件を満たさない解への遷移が生じた場合のコスト項の変化量の期待値を計算し、期待値に基づいて制約項の係数を決定し、決定した係数を用いてエネルギー関数の情報を出力する。
【選択図】図１
【特許請求の範囲】
【請求項1】
制約条件を含む組合せ最適化問題を、コスト項と前記制約条件に対応する制約項との和で表されるイジングモデルのエネルギー関数に基づいて計算する情報処理システムであって、
前記組合せ最適化問題の問題情報を取得し、前記問題情報に基づいて、前記制約条件を満たす解から前記制約条件を満たさない解への遷移が生じた場合の前記コスト項の変化量の期待値を計算し、前記期待値に基づいて前記制約項の係数を決定し、決定した前記係数を用いて前記エネルギー関数の情報を出力する処理部と、
出力された前記情報に基づいて、前記イジングモデルの基底状態を探索する探索部と、
を有する情報処理システム。
前記組合せ最適化問題の問題情報を取得し、前記問題情報に基づいて、前記制約条件を満たす解から前記制約条件を満たさない解への遷移が生じた場合の前記コスト項の変化量の期待値を計算し、前記期待値に基づいて前記制約項の係数を決定し、決定した前記係数を用いて前記エネルギー関数の情報を出力する処理部と、
出力された前記情報に基づいて、前記イジングモデルの基底状態を探索する探索部と、
を有する情報処理システム。
【請求項2】
前記組合せ最適化問題は二次割当て問題であり、
前記処理部は、前記二次割当て問題の前記問題情報を表す第１の行列と第２の行列とのクロネッカー積によって得られる第３の行列を計算し、前記第３の行列に含まれる要素の第１の平均値に基づいて前記期待値を計算する、請求項１に記載の情報処理システム。
前記処理部は、前記二次割当て問題の前記問題情報を表す第１の行列と第２の行列とのクロネッカー積によって得られる第３の行列を計算し、前記第３の行列に含まれる要素の第１の平均値に基づいて前記期待値を計算する、請求項１に記載の情報処理システム。
【請求項3】
前記処理部は、前記期待値に調整係数を乗ずることで、前記係数を計算する、請求項１に記載の情報処理システム。
【請求項4】
前記処理部が、前記期待値に前記調整係数を乗じて得られた前記係数を用いて定式化した前記エネルギー関数の前記情報を出力し、
前記探索部が、前記情報に基づいて行った前記基底状態の探索結果を出力し、
前記処理部が、前記調整係数を増加または減少させる、
処理を前記調整係数が上限または下限に達するまで繰返し、
前記処理部は、前記調整係数が上限または前記下限に達した場合、前記探索結果に基づいて、前記調整係数を決定する、
請求項３に記載の情報処理システム。
前記探索部が、前記情報に基づいて行った前記基底状態の探索結果を出力し、
前記処理部が、前記調整係数を増加または減少させる、
処理を前記調整係数が上限または下限に達するまで繰返し、
前記処理部は、前記調整係数が上限または前記下限に達した場合、前記探索結果に基づいて、前記調整係数を決定する、
請求項３に記載の情報処理システム。
【請求項5】
前記第３の行列は、前記第１の行列または前記第２の行列のうちの一方の行列に含まれる各要素と、他方の行列との積による複数の小行列を含み、
前記処理部は、前記複数の小行列のそれぞれにおいて、前記複数の小行列のそれぞれに含まれる要素の第２の平均値と前記第１の平均値との差分値を、前記複数の小行列のそれぞれに含まれる前記要素の値に加算または減算することで、前記第１の平均値と前記第２の平均値とを一致させる、
請求項２に記載の情報処理システム。
前記処理部は、前記複数の小行列のそれぞれにおいて、前記複数の小行列のそれぞれに含まれる要素の第２の平均値と前記第１の平均値との差分値を、前記複数の小行列のそれぞれに含まれる前記要素の値に加算または減算することで、前記第１の平均値と前記第２の平均値とを一致させる、
請求項２に記載の情報処理システム。
【請求項6】
制約条件を含む組合せ最適化問題を、コスト項と前記制約条件に対応する制約項との和で表されるイジングモデルのエネルギー関数に基づいて計算する情報処理方法であって、
処理部が、前記組合せ最適化問題の問題情報を取得し、前記問題情報に基づいて、前記制約条件を満たす解から前記制約条件を満たさない解への遷移が生じた場合の前記コスト項の変化量の期待値を計算し、前記期待値に基づいて前記制約項の係数を決定し、決定した前記係数を用いて前記エネルギー関数の情報を出力し、
探索部が、出力された前記情報に基づいて、前記イジングモデルの基底状態を探索する、
情報処理方法。
処理部が、前記組合せ最適化問題の問題情報を取得し、前記問題情報に基づいて、前記制約条件を満たす解から前記制約条件を満たさない解への遷移が生じた場合の前記コスト項の変化量の期待値を計算し、前記期待値に基づいて前記制約項の係数を決定し、決定した前記係数を用いて前記エネルギー関数の情報を出力し、
探索部が、出力された前記情報に基づいて、前記イジングモデルの基底状態を探索する、
情報処理方法。
【請求項7】
制約条件を含む組合せ最適化問題を、コスト項と前記制約条件に対応する制約項との和で表されるイジングモデルのエネルギー関数に基づいて計算する情報処理システムに用いられるコンピュータに、
前記組合せ最適化問題の問題情報を取得し、前記問題情報に基づいて、前記制約条件を満たす解から、前記制約条件を満たさない解への遷移が生じた場合の前記コスト項の変化量の期待値を計算し、
前記期待値に基づいて前記制約項の係数を決定し、
決定した前記係数を用いて前記エネルギー関数の情報を、前記情報に基づいて前記イジングモデルの基底状態を探索する探索部に出力する、
処理を実行させるプログラム。
前記組合せ最適化問題の問題情報を取得し、前記問題情報に基づいて、前記制約条件を満たす解から、前記制約条件を満たさない解への遷移が生じた場合の前記コスト項の変化量の期待値を計算し、
前記期待値に基づいて前記制約項の係数を決定し、
決定した前記係数を用いて前記エネルギー関数の情報を、前記情報に基づいて前記イジングモデルの基底状態を探索する探索部に出力する、
処理を実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、情報処理システム、情報処理方法及びプログラムに関する。
【背景技術】
【0002】
ノイマン型コンピュータが不得意とする多変数の組合せ最適化問題を、磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えて計算するイジングマシン（ボルツマンマシンとも呼ばれる）がある。イジングモデルに置き換えられた問題を実用的な時間で解く手法には、たとえば、シミュレーテッドアニーリング（ＳＡ：Simulated Annealing）法やレプリカ交換法などのマルコフ連鎖モンテカルロ（ＭＣＭＣ：Markov Chain Monte Carlo）法がある。組合せ最適化問題は、複数の状態変数を含むエネルギー関数により定式化される。たとえば、情報処理装置は、ＭＣＭＣ法を用いて状態変数の値を変化させることによる状態遷移を繰り返し試行することで、エネルギー関数の値が最小となるイジングモデルの基底状態を探索する。基底状態は組合せ最適化問題の最適解に対応する。
【0003】
エネルギー関数は、コスト項と制約項との和で表すことができる（たとえば、特許文献１，２参照）。エネルギー関数（Ｅ（ｘ））は、たとえば、以下の式（１）で表せる。
【0004】
【数1】
【0005】
式（１）において、ｘはｘ＝（ｘ１，ｘ２，…，ｘｎ）のバイナリ変数である状態変数によるベクトルである。Ｃ（ｘ）はコスト項であり、最小化したい目的関数である。αＰ（ｘ）は制約項である。組合せ最適化問題が満たすべき制約条件を満たすとき（制約充足時）には、Ｐ（ｘ）＝０となり、制約条件を満たさないとき（制約違反時）には、Ｐ（ｘ）＞０となる。αは、制約係数であり、０より大きい値をもつ。このような制約項をコスト項に加えることで、制約条件を満たさない解（制約違反解）が選ばれにくくなる。
【0006】
なお、エネルギー関数の符号を変えれば、エネルギー関数の値が最大となるイジングモデルの状態を探索することもできる。
組合せ最適化問題の実用的な問題の例として、二次割当問題（ＱＡＰ：Quadratic Assignment Problem）が挙げられる。ＱＡＰは、ｎ個の割当要素（施設など）をｎ個の割当先（場所など）に割り当てる際に、割当要素間のフロー量（施設間での物資の輸送量などのコスト）と、各割当要素が割り当てられる割当先間の距離との積の総和を最小化するような割当を求める問題である。割当要素間のフロー量と、各割当先間の距離とはそれぞれ行列を用いて表せる。また、ＱＡＰには、全ての割当要素を、互いに異なる割当先に割り当てる、という制約条件が課される。
組合せ最適化問題の実用的な問題の例として、二次割当問題（ＱＡＰ：Quadratic Assignment Problem）が挙げられる。ＱＡＰは、ｎ個の割当要素（施設など）をｎ個の割当先（場所など）に割り当てる際に、割当要素間のフロー量（施設間での物資の輸送量などのコスト）と、各割当要素が割り当てられる割当先間の距離との積の総和を最小化するような割当を求める問題である。割当要素間のフロー量と、各割当先間の距離とはそれぞれ行列を用いて表せる。また、ＱＡＰには、全ての割当要素を、互いに異なる割当先に割り当てる、という制約条件が課される。
【0007】
ところで、上記のような制約係数（式（１）の場合はα）は、小さすぎると制約違反解が発生しやすくなり、逆に大きすぎると制約充足解間のエネルギー障壁が大きくなって解空間の探索に支障をきたす。また、適切な制約係数の値は個々の問題に依存して大きく異なる。従来、イジングマシンによる組合せ最適化問題の解の探索結果に基づいて、制約係数を調整する処理を繰り返すことで、制約係数を最適化していく手法が提案されている。
【先行技術文献】
【特許文献】
【0008】
【特許文献1】国際公開第２０２１－０４４５１６号
【特許文献2】特開２００４－１１０８３１号公報
【発明の概要】
【発明が解決しようとする課題】
【0009】
しかし、従来の手法で、高品質の解を得るために広範囲での探索ができるような制約係数を得るには多数回の調整が必要となる場合があり、制約係数の決定に時間がかかる。
１つの側面では、本発明は、制約係数の決定にかかる時間を短縮可能な情報処理システム、情報処理方法及びプログラムを提供することを目的とする。
１つの側面では、本発明は、制約係数の決定にかかる時間を短縮可能な情報処理システム、情報処理方法及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【0010】
１つの実施態様では、制約条件を含む組合せ最適化問題を、コスト項と前記制約条件に対応する制約項との和で表されるイジングモデルのエネルギー関数に基づいて計算する情報処理システムであって、前記組合せ最適化問題の問題情報を取得し、前記問題情報に基づいて、前記制約条件を満たす解から前記制約条件を満たさない解への遷移が生じた場合の前記コスト項の変化量の期待値を計算し、前記期待値に基づいて前記制約項の係数を決定し、決定した前記係数を用いて前記エネルギー関数の情報を出力する処理部と、出力された前記情報に基づいて、前記イジングモデルの基底状態を探索する探索部と、を有する情報処理システムが提供される。
【0011】
また、１つの実施態様では、情報処理方法が提供される
また、１つの実施態様では、プログラムが提供される。
また、１つの実施態様では、プログラムが提供される。
【発明の効果】
【0012】
１つの側面では、本発明は、制約係数の決定にかかる時間を短縮できる。
【図面の簡単な説明】
【0013】
【図1】第１の実施の形態の情報処理システムの一例を示す図である。
【図2】第２の実施の形態の情報処理システムのハードウェア例を示す図である。
【図3】情報処理システムの機能例を示す図である。
【図4】ＱＡＰにおける制約違反の有無と制約違反エネルギーの例を示す図である。
【図5】状態とエネルギー関数の値との関係の例を示す図である。
【図6】情報処理システムの一例の処理の流れを示すフローチャートである。
【図7】行列展開例を示す図である。
【図8】ある割当状態のときにＱｉｊから選択される要素ｑｉｊの変化の例を示す図である。
【図9】比較例の情報処理システムの処理の流れを示すフローチャートである。
【図10】比較例における制約係数の値とインスタンスの計算結果との関係を示す図である。
【図11】本実施の形態における調整係数の値とインスタンスの計算結果との関係を示す図である。
【図12】調整係数の決定処理の一例の流れを示すフローチャートである。
【図13】調整係数の値と問題サイズの異なる同種のインスタンスの計算結果との関係を示す図である。
【図14】式（１６）における２次項に足し合わされるｘｉｋｘｊｋと、１次項とを行列で表した図である。
【図15】式（１７）における２次項に足し合わされるｘｉｋｘｉｌと、１次項とを行列で表した図である。
【図16】制約項のＱＵＢＯ行列の例を示す図である。
【図17】エネルギー関数全体のＱＵＢＯ行列の例を示す図である。
【図18】第３の実施の形態の情報処理システムの一例の処理の流れを示すフローチャートである。
【図19】小行列ｄｋｌＦの要素の小行列ｆｉｊＤへの配分例を示す図である。
【図20】行列変調の一例の流れを示すフローチャートである。
【図21】行列変調を用いた場合の調整係数の値とインスタンスの計算結果との関係を示す図である。
【発明を実施するための形態】
【0014】
以下、発明を実施するための形態を、図面を参照しつつ説明する。
（制約係数の決定方針）
イジングモデルのエネルギー関数は、前述の式（１）で表せる。ここで、制約係数であるαは、制約違反を防止するため、制約違反解が制約充足解よりも大きくなるような値とすればよい。しかし、制約充足解は、Ｃ（ｘ）の値が最小である最適解から、制約条件が充足されていても、Ｃ（ｘ）の値が大きなものまで幅があり、どの値を基準にするのが適当かはわからない。また、そもそも最適解の値が明らかではない。前述のようにイジングマシンによる組合せ最適化問題の解の探索結果に基づいて、制約係数を調整する処理を繰り返すことで、制約係数を最適化していくことが提案されているが、調整に時間がかかる。そこで、本実施の形態では、制約充足解から制約違反解への遷移時のＣ（ｘ）の値の変化量に着目して、αの決定が行われる。
（制約係数の決定方針）
イジングモデルのエネルギー関数は、前述の式（１）で表せる。ここで、制約係数であるαは、制約違反を防止するため、制約違反解が制約充足解よりも大きくなるような値とすればよい。しかし、制約充足解は、Ｃ（ｘ）の値が最小である最適解から、制約条件が充足されていても、Ｃ（ｘ）の値が大きなものまで幅があり、どの値を基準にするのが適当かはわからない。また、そもそも最適解の値が明らかではない。前述のようにイジングマシンによる組合せ最適化問題の解の探索結果に基づいて、制約係数を調整する処理を繰り返すことで、制約係数を最適化していくことが提案されているが、調整に時間がかかる。そこで、本実施の形態では、制約充足解から制約違反解への遷移時のＣ（ｘ）の値の変化量に着目して、αの決定が行われる。
【0015】
なお、温度パラメータを用いるＭＣＭＣ法では、温度パラメータの値の設定にもよるが、Ｃ（ｘ）の値が増加するような制約充足解から制約違反解への遷移が受け入れられる確率は低い。このため、以下では、制約充足解から制約違反解への遷移によりＣ（ｘ）の値が減少するケースを考慮する。
【0016】
制約充足解から制約違反解への遷移によりＣ（ｘ）の値が減少する場合、制約違反により加算されるαＰ（ｘ）が、Ｃ（ｘ）の変化量（減少量）よりも小さければその遷移により式（１）のＥ（ｘ）が小さくなる。このため、高確率でその制約違反解への遷移が起こりうる。したがって、制約違反を抑止するためにはαＰ（ｘ）は、少なくともＣ（ｘ）の減少量を補償する程度の値であればよい。次に、Ｃ（ｘ）の減少量の期待値を計算する方法を説明する。
【0017】
（Ｃ（ｘ）の減少量の期待値の計算方法）
式（１）において、Ｃ（ｘ）を、ＱＵＢＯ（Quadratic Unconstrained Binary Optimization）で表すと、以下の式（２）のようになる。
式（１）において、Ｃ（ｘ）を、ＱＵＢＯ（Quadratic Unconstrained Binary Optimization）で表すと、以下の式（２）のようになる。
【0018】
【数2】
【0019】
式（２）において、Ｑｉｊは結合係数行列、ｈｉはバイアス係数ベクトル、ｃは定数項である。
Ｑｉｊの次数をｎとするとＱｉｊの要素ｑｉｊ（ｉ＜ｊ）の数は、ｎ２個であるが、ｑｉｊ＝ｑｊｉ、ｑｉｉ＝０（つまり対角要素は０）であるため、有効なｑｉｊの数は、Ｎ＝ｎ（ｎ－１）／２である。以下、上記有効なＮ個のｑｉｊを、結合係数行列であるＱｉｊの上三角の要素（対角要素を除く）とし、ｑｓ（ｓ＝１，…，Ｎ）と表記する。ｑｓの総和は、以下の式（３）で表せ、ｑｓの平均値（μｑ）は、以下の式（４）で表せる。
Ｑｉｊの次数をｎとするとＱｉｊの要素ｑｉｊ（ｉ＜ｊ）の数は、ｎ２個であるが、ｑｉｊ＝ｑｊｉ、ｑｉｉ＝０（つまり対角要素は０）であるため、有効なｑｉｊの数は、Ｎ＝ｎ（ｎ－１）／２である。以下、上記有効なＮ個のｑｉｊを、結合係数行列であるＱｉｊの上三角の要素（対角要素を除く）とし、ｑｓ（ｓ＝１，…，Ｎ）と表記する。ｑｓの総和は、以下の式（３）で表せ、ｑｓの平均値（μｑ）は、以下の式（４）で表せる。
【0020】
【数3】
【0021】
【数4】
【0022】
また、ｈｉの要素数はｎであり、ｈｉの平均値（μｈ）は、以下の式（５）で表せる。
【0023】
【数5】
【0024】
式（２）のＱｉｊｘｉｘｊは、ｘｉｘｊ＝１の場合に０以外の値をもつ。つまり、ｘｉｘｊ＝１の場合にＱｉｊのｉ列とｊ行の交差する要素が選択される。ある解ｘに含まれるｘｉ（ｉ＝１，…，ｎ）において、ｘｉ＝１である状態変数の個数をｐ個とすると、選択されるＱｉｊの要素の数は、ｐ（ｐ－１）／２である。
【0025】
式（２）のｈｉｘｉは、ｘｉ＝１の場合に０以外の値をもつ。つまり、ｘｉ＝１の場合、ｈｉの要素が選択される。ｘｉ＝１である状態変数の個数をｐ個とすると、選択されるｈｉの要素の数は、ｐ個である。
【0026】
値が１であるｐ個のｘｉのうちの１つが、１から０に反転する（非選択になる）とき、Ｑｉｊのうち非選択になる要素の数はｐ－１個、ｈｉのうち非選択になる要素の数は１個である。このＱｉｊとｈｉの非選択になる要素の和がＣ（ｘ）の減少量となる。
【0027】
上記のｑｓ（ｓ＝１，…，Ｎ）から要素を抽出する試行を考えると、抽出される要素の値Ｘ＝ｑ１，…，ｑＮは確率変数である。ただし、上述の非選択になるｐ－１個の要素の抽出は非復元抽出であり、また、非選択になるｐ－１個の要素はＱｉｊの同じ行または列に存在し、さらに反転するｘｉの選択も無作為ではなくｘの関数となっている。したがって、確率変数列｛Ｘ１，…，Ｘｐ－１｝は独立同分布にしたがうわけではない。
【0028】
しかし、非復元抽出はｎが大きい場合には近似的に復元抽出とみなせる。また、反転するｘｉの選択においては、イジングマシンによる求解は確率的過程を含む。そこで、近似としてＸを独立同分布（ｉｉｄ）にしたがう確率変数と仮定すると、非選択になるｐ－１個の要素で構成される確率変数列｛Ｘ１，…，Ｘｐ－１｝の和Ｓの期待値は、μｑ（１→０）＝（ｐ－１）μｑとなる。
【0029】
同様にバイアス係数の集合｛ｈｓ（ｓ＝１，…，Ｎ）｝から要素を抽出する試行を考える。抽出される要素の値Ｙ＝ｈ１，…，ｈｐをｉｉｄにしたがう確率変数と仮定すると、非選択になる要素で構成される確率変数列｛Ｙ１｝（上記のように要素数は１個）の和の期待値は、μｈ（１→０）＝μｈとなる。
【0030】
以上より、Ｃ（ｘ）の減少量の期待値は、Δ－Ｃ（ｘ）＝μｑ（１→０）＋μｈ（１→０）＝（ｐ－１）μｑ＋μｈとなる。このようなΔ－Ｃ（ｘ）を制約係数αとして用いれば、制約充足解から制約違反解への遷移によりＣ（ｘ）の値が減少する場合、制約違反により加算されるαＰ（ｘ）を、Ｃ（ｘ）の減少量を補償する程度の値とすることができる。
【0031】
（ＱＡＰへの適用）
ＱＡＰの目的関数ではｈｉとｃは恒等的に０であるため、Ｃ（ｘ）は以下の式（６）で表せる。
ＱＡＰの目的関数ではｈｉとｃは恒等的に０であるため、Ｃ（ｘ）は以下の式（６）で表せる。
【0032】
【数6】
【0033】
ＱＡＰでは、割当要素間のフロー量と、各割当先間の距離とはそれぞれ行列を用いて表せる（後述の図１参照）。ｎ個の割当要素をｎ個の割当先に割り当てる場合のＣ（ｘ）は、以下の式（７）で表せる。
【0034】
【数7】
【0035】
式（７）において、ｆｉｊは、割当要素間のフロー量を表すフロー行列Ｆのｉ列ｊ行のフロー量（ｉ番目の割当要素とｊ番目の割当要素との間のフロー量）を表す。ｄｋｌは、割当先間の距離を表す距離行列Ｄのｋ列ｌ行の距離（ｋ番目の割当先とｌ番目の割当先との間の距離）を表す。ｘｉｋは、割当状態を表す割当行列ｘ（後述の図１参照）のｉ列ｋ行の割当を表す。ｉ番目の割当要素がｋ番目の割当先に割り当てられている場合には、ｘｉｋ＝１、ｉ番目の割当要素がｋ番目の割当先に割り当てられていない場合には、ｘｉｋ＝０となる。
【0036】
なお、式（１）のＰ（ｘ）は、ＱＡＰでは式（８）のように表すことができる。
【0037】
【数8】
【0038】
右辺の１項目は、同じｋでｘｉｋ＝１となるｉを１個にするための制約項であり、１つの割当先には１つの割当要素が割り当てられるという制約を表す。右辺の２項目は、同じｉでｘｉｋ＝１となるｋを１個にするための制約項であり、１つの割当要素は１つの割当先に割り当てられるという制約を表す。α≠βとすることもあり得るが、本質的にＱＡＰではフロー行列Ｆと距離行列Ｄは対称行列であるため、β＝αとし、まとめてαＰ（ｘ）として制約項を表すこともできる。
【0039】
式（６）のＱｉｊは、フロー行列Ｆと距離行列Ｄのクロネッカー積によって、以下の式（９）のような行列Ｑで表せる。行列Ｑは、ｎ×ｎ個の小行列ｆｉｊＤからなる。
【0040】
【数9】
【0041】
式（９）を展開すると、式（１０）に示すような次数がｎ２の行列Ｑが得られる。
【0042】
【数10】
【0043】
また、割当行列ｘは、以下の式（１１）のように、要素数がｎ２個の列ベクトルで表せる。
【0044】
【数11】
【0045】
式（１０）のように、行列Ｑの要素ｑｉｊ（ｉ＜ｊ）の数は、ｎ４個であるが、ｑｉｊ＝ｑｊｉ、ｑｉｉ＝０（つまり対角要素は０）であるため、有効なｑｉｊの数は、Ｎ＝ｎ２（ｎ２－１）／２である。有効なＮ個のｑｉｊを、行列Ｑの上三角の要素（対角要素を除く）とし、ｑｓ（ｓ＝１，…，Ｎ）と表記すると、ｑｓの総和は、前述の式（３）で表せ、ｑｓの平均値は、前述の式（４）で表せる。
【0046】
Ｃ（ｘ）の減少量の期待値は、前述のようにΔ－Ｃ（ｘ）＝μｑ（１→０）＋μｈ（１→０）＝（ｐ－１）μｑ＋μｈと表せるが、ＱＡＰの場合、μｈ≡０である。また、ＱＡＰの場合、制約充足解で１となるｘｉの個数ｐは、ｐ＝ｎである。したがって、Δ－Ｃ（ｘ）は、Δ－Ｃ（ｘ）＝（ｐ－１）μｑ＝（ｎ－１）μｑ＝（ｎ－１）Ｓ／Ｎ＝２Ｓ／ｎ２（ｎ＋１）となる。
【0047】
このようなΔ－Ｃ（ｘ）を制約係数αとして用いれば、制約充足解から制約違反解への遷移によりＣ（ｘ）の値が減少する場合、制約違反により加算されるαＰ（ｘ）を、Ｃ（ｘ）の減少量を補償する程度の値とすることができる。
【0048】
（Δ－Ｃ（ｘ）の他の例）
式（９）に示したように、行列Ｑは、ｎ×ｎ個の小行列ｆｉｊＤからなる。行列Ｑの上三角で対角要素（ｆｉｊＤ（ｉ＝ｊ））を含まない小行列の数は、ｎ（ｎ－１）／２個であり、各小行列において対角要素（ｆｉｊｄｋｌ（ｋ＝ｌ））を含まない要素数はｎ（ｎ－１）個である。この場合、有効なｑｉｊ（ｉ＜ｊ）の数は、Ｎ’＝ｎ（ｎ－１）×ｎ（ｎ－１）／２＝ｎ２（ｎ－１）２／２である。上記有効なＮ’個のｑｉｊを、ｑｓ（ｓ＝１，…，Ｎ’）と表記すると、ｑｓの総和Ｓ’は、以下の式（１２）で表せる。
式（９）に示したように、行列Ｑは、ｎ×ｎ個の小行列ｆｉｊＤからなる。行列Ｑの上三角で対角要素（ｆｉｊＤ（ｉ＝ｊ））を含まない小行列の数は、ｎ（ｎ－１）／２個であり、各小行列において対角要素（ｆｉｊｄｋｌ（ｋ＝ｌ））を含まない要素数はｎ（ｎ－１）個である。この場合、有効なｑｉｊ（ｉ＜ｊ）の数は、Ｎ’＝ｎ（ｎ－１）×ｎ（ｎ－１）／２＝ｎ２（ｎ－１）２／２である。上記有効なＮ’個のｑｉｊを、ｑｓ（ｓ＝１，…，Ｎ’）と表記すると、ｑｓの総和Ｓ’は、以下の式（１２）で表せる。
【0049】
【数12】
【0050】
ここで、ＱＡＰでは、行列Ｑの無効な要素の値は０であるため、Ｓ’＝Ｓである。このためｑｓの平均値は、以下の式（１３）で表せる。
【0051】
【数13】
【0052】
したがって、Δ－Ｃ（ｘ）は、Δ－Ｃ（ｘ）＝（ｎ－１）Ｓ／Ｎ’＝２（ｎ－１）Ｓ／ｎ２（ｎ－１）２＝２Ｓ／ｎ２（ｎ－１）となる。
このようなΔ－Ｃ（ｘ）を制約係数αとして用いても、制約充足解から制約違反解への遷移によりＣ（ｘ）の値が減少する場合、制約違反により加算されるαＰ（ｘ）を、Ｃ（ｘ）の減少量を補償する程度の値とすることができる。
このようなΔ－Ｃ（ｘ）を制約係数αとして用いても、制約充足解から制約違反解への遷移によりＣ（ｘ）の値が減少する場合、制約違反により加算されるαＰ（ｘ）を、Ｃ（ｘ）の減少量を補償する程度の値とすることができる。
【0053】
（第１の実施の形態）
図１は、第１の実施の形態の情報処理システムの一例を示す図である。
情報処理システム１０は、組合せ最適化問題を、式（１）に示したようなイジングモデルのエネルギー関数に基づいて計算する。
図１は、第１の実施の形態の情報処理システムの一例を示す図である。
情報処理システム１０は、組合せ最適化問題を、式（１）に示したようなイジングモデルのエネルギー関数に基づいて計算する。
【0054】
情報処理システム１０は、記憶部１１、処理部１２及び探索部１３を有する。
記憶部１１は、ＤＲＡＭ（Dynamic Random Access Memory）などの揮発性記憶装置でもよいし、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性記憶装置でもよい。記憶部１１は、処理部１２の処理に用いられる各種のデータを記憶する。また、記憶部１１は、組合せ最適化問題の問題情報を記憶する。たとえば、最適化問題がＱＡＰの場合、問題情報は、各割当要素間のフロー量と、各割当先間の距離とを含む。
記憶部１１は、ＤＲＡＭ（Dynamic Random Access Memory）などの揮発性記憶装置でもよいし、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性記憶装置でもよい。記憶部１１は、処理部１２の処理に用いられる各種のデータを記憶する。また、記憶部１１は、組合せ最適化問題の問題情報を記憶する。たとえば、最適化問題がＱＡＰの場合、問題情報は、各割当要素間のフロー量と、各割当先間の距離とを含む。
【0055】
処理部１２は、ＣＰＵ（Central Processing Unit）、ＤＳＰ（Digital Signal Processor）、ＡＳＩＣ（Application Specific Integrated Circuit）、ＦＰＧＡ（Field Programmable Gate Array）などを含み得る。処理部１２はプログラムを実行するプロセッサでもよい。「プロセッサ」には、複数のプロセッサの集合（マルチプロセッサ）が含まれ得る。
【0056】
処理部１２は、組合せ最適化問題の問題情報を取得し、問題情報に基づいて、組合せ最適化問題をエネルギー関数で定式化し、エネルギー関数の情報を探索部１３に出力する。処理部１２は、探索部１３で求められた解を取得し、取得した解を外部に出力する。たとえば、処理部１２は、情報処理システム１０に接続された表示装置に解の内容を表示させてもよいし、ネットワークを介して他の情報処理装置に解を送信してもよい。
【0057】
探索部１３は、処理部１２からエネルギー関数の情報を受け付け、エネルギー関数の情報に基づいて、イジングモデルの基底状態、すなわち、エネルギー関数の値を最小化する解を探索する。探索部１３は、探索により得られた解を処理部１２に出力する。探索部１３は、専用のハードウェアにより実現されうる。たとえば、ＦＰＧＡ、ＡＳＩＣ、ＧＰＵ（Graphics Processing Unit）などの集積回路を用いて実現される探索回路が、探索部１３として機能してもよい。ただし、探索部１３は、プロセッサが所定のプログラムを実行することで実現されてもよい。探索部１３が用いる探索手法には、たとえば、ＳＡ法、レプリカ交換法及びＳＱＡ（Simulated Quantum Annealing）などがある。
【0058】
以下、第１の実施の形態の情報処理システム１０の動作例を説明する。
処理部１２は、組合せ最適化問題をエネルギー関数で定式化するにあたり、制約項の係数（制約係数）であるαの値を決定する。
処理部１２は、組合せ最適化問題をエネルギー関数で定式化するにあたり、制約項の係数（制約係数）であるαの値を決定する。
【0059】
まず、処理部１２は、記憶部１１に記憶された組合せ最適化問題の問題情報を取得する（読み出す）。図１には、ＱＡＰの問題情報の例が示されている。
ＱＡＰの問題情報は、各割当要素間のフロー量と、各割当先間の距離とを含む。各割当要素間のフロー量は、ｎ行ｎ列のフロー行列Ｆで表され、各割当先間の距離は、ｎ行ｎ列の距離行列Ｄで表される。
ＱＡＰの問題情報は、各割当要素間のフロー量と、各割当先間の距離とを含む。各割当要素間のフロー量は、ｎ行ｎ列のフロー行列Ｆで表され、各割当先間の距離は、ｎ行ｎ列の距離行列Ｄで表される。
【0060】
なお、図１には、割当行列ｘの例も示されている。ｘｉｋは、割当状態を表す割当行列ｘのｉ列ｋ行の割当を表す。ｉ番目の割当要素がｋ番目の割当先に割り当てられている場合には、ｘｉｋ＝１、ｉ番目の割当要素がｋ番目の割当先に割り当てられていない場合には、ｘｉｋ＝０となる。なお、割当行列ｘは、前述の式（１１）のように、要素数がｎ２個の列ベクトルで表すこともできる。
【0061】
処理部１２は、取得した問題情報に基づいて、組合せ最適化問題の制約条件を満たす解（制約充足解）から、制約条件を満たさない解（制約違反解）への遷移が生じた場合のコスト項の変化量の期待値を計算する。
【0062】
コスト項が式（２）のように表される場合、前述の（Ｃ（ｘ）の減少量の期待値の計算方法）のように、変化量（減少量）の期待値は、Δ－Ｃ（ｘ）＝（ｐ－１）μｑ＋μｈと表される式により計算できる。なお、前述のように、ｐは、ある解に含まれるｎ個の状態変数ｘｉ（ｉ＝１，…，ｎ）において、ｘｉ＝１である状態変数の個数である。μｑは、結合係数行列であるＱｉｊにおいて有効なＮ個の要素ｑｓ（ｓ＝１，…，Ｎ）の平均値であり、μｈは、ｈｉ（バイアス係数ベクトル）のｎ個の要素の平均値である。
【0063】
記憶部１１には、制約充足解が初期解として記憶されており、その初期解に基づいて、処理部１２がｐの値を決定してもよいし、ｐの値自体が予め記憶部１１に記憶されていてもよい。μｑやμｈは問題情報に基づいて、前述の式（３）～式（５）により計算される。組合せ最適化問題がＴＳＰ（Traveling Salesman Problem）などの場合、上記のΔ－Ｃ（ｘ）＝（ｐ－１）μｑ＋μｈにより、Δ－Ｃ（ｘ）の計算が可能である。
【0064】
図１には、組合せ最適化問題がＱＡＰの場合のΔ－Ｃ（ｘ）の計算例が示されている。ＱＡＰのコスト項Ｃ（ｘ）は、前述の式（６）のように表せる。処理部１２は、前述の式（９）に示したようにフロー行列Ｆと距離行列Ｄのクロネッカー積によって行列Ｑを計算する。そして、処理部１２は、式（１０）に示したような行列Ｑに含まれる各要素の平均値μｑに基づいて、Δ－Ｃ（ｘ）を計算する。なお、式（６）、式（９）、式（１０）は、図１にも示されている。
【0065】
ＱＡＰの場合、前述のように、行列Ｑの上三角の要素（対角要素を除く）である有効な要素の平均値を求めればよい。有効な要素の数は、Ｎ＝ｎ２（ｎ２－１）／２であるから、有効な要素の値の和ＳをＮで割った、Ｓ／Ｎ＝２Ｓ／ｎ２（ｎ２－１）が平均値μｑである。さらに、組合せ最適化問題がＱＡＰの場合、ｐ＝ｎであるため、Δ－Ｃ（ｘ）＝２Ｓ／ｎ２（ｎ＋１）となる。
【0066】
なお、処理部１２は、前述の（Δ－Ｃ（ｘ）の他の例）で示したように、Δ－Ｃ（ｘ）＝２Ｓ／ｎ２（ｎ－１）と表される式によりΔ－Ｃ（ｘ）を計算してもよい。
また、ＱＡＰの問題サイズ（ｎ）が大きくなると、Δ－Ｃ（ｘ）の計算式に含まれるｎ＋１やｎ－１はほぼ等しくなり、ｎともほぼ等しくなる。このため、ｎ＋１やｎ－１の代わりに、ｎを用いてもΔ－Ｃ（ｘ）の値はほぼ同様の値となる。したがって、上記のΔ－Ｃ（ｘ）の計算式は、問題サイズに応じて適宜変更可能である。
また、ＱＡＰの問題サイズ（ｎ）が大きくなると、Δ－Ｃ（ｘ）の計算式に含まれるｎ＋１やｎ－１はほぼ等しくなり、ｎともほぼ等しくなる。このため、ｎ＋１やｎ－１の代わりに、ｎを用いてもΔ－Ｃ（ｘ）の値はほぼ同様の値となる。したがって、上記のΔ－Ｃ（ｘ）の計算式は、問題サイズに応じて適宜変更可能である。
【0067】
処理部１２は、計算したΔ－Ｃ（ｘ）に基づいて制約係数αを決定する。処理部１２は、たとえば、α＝Δ－Ｃ（ｘ）と決定してもよいし、Δ－Ｃ（ｘ）に、後述の調整係数を乗じた値をαとして決定してもよい。
【0068】
その後、処理部１２は、決定したαを用いて組合せ最適化問題を定式化したエネルギー関数の情報を出力する。出力されるエネルギー関数の情報は、たとえば、式（１）のエネルギー関数をＱＵＢＯ形式で表した行列データ（ＱＵＢＯ行列とも呼ばれる）である。
【0069】
探索部１３による探索は、たとえば、次のように実行される。一例として、ＳＡ法またはレプリカ交換法を用いる場合を考える。
探索部１３は、複数の状態変数のうちの１つの状態変数の値が変化することによるエネルギー関数の値の変化量を、複数の状態変数の各々について計算し、エネルギー関数の値が小さくなる変化が優先される形で確率的に受け入れる。このとき、最急降下法では、局所解に陥った場合に脱出できなくなる。そこで、探索部１３は、ある状態変数を変化させることによる、イジングモデルのある状態から次の状態への遷移確率の決定に、メトロポリス法やギブス法を用いる。すなわち、探索部１３は、エネルギー関数の値が大きくなる変化についても、エネルギー関数の値の変化量とノイズ値との比較に応じて確率的に許容する。ノイズ値は、温度パラメータの値や乱数に基づいて求められる。温度パラメータの値が大きい程、ノイズ値の振幅が大きくなる。ノイズ値の振幅が大きい程、エネルギー関数の値の増加量が大きい状態遷移が許容されやすくなる。
探索部１３は、複数の状態変数のうちの１つの状態変数の値が変化することによるエネルギー関数の値の変化量を、複数の状態変数の各々について計算し、エネルギー関数の値が小さくなる変化が優先される形で確率的に受け入れる。このとき、最急降下法では、局所解に陥った場合に脱出できなくなる。そこで、探索部１３は、ある状態変数を変化させることによる、イジングモデルのある状態から次の状態への遷移確率の決定に、メトロポリス法やギブス法を用いる。すなわち、探索部１３は、エネルギー関数の値が大きくなる変化についても、エネルギー関数の値の変化量とノイズ値との比較に応じて確率的に許容する。ノイズ値は、温度パラメータの値や乱数に基づいて求められる。温度パラメータの値が大きい程、ノイズ値の振幅が大きくなる。ノイズ値の振幅が大きい程、エネルギー関数の値の増加量が大きい状態遷移が許容されやすくなる。
【0070】
探索部１３は、このような処理に基づく状態遷移を、探索開始から探索終了条件が満たされるまで繰り返し試行し、探索終了条件が満たされたときの複数の状態変数の値の組を解とする。たとえば、探索終了条件は、所定の単位時間が経過したか、あるいは、状態変数の変化の試行回数が所定の単位回数に到達したかなどにより判定される。
【0071】
以上のように、情報処理システム１０の処理部１２は、制約充足解から制約違反解への遷移が生じた場合のＣ（ｘ）の変化量の期待値Δ－Ｃ（ｘ）に基づいて、制約係数αを決定する。これにより、制約充足解から制約違反解への遷移によりＣ（ｘ）の値が減少する場合、制約違反により加算されるαＰ（ｘ）を、Ｃ（ｘ）の減少量を補償する程度の値とすることができる。したがって、適切なαを得るために多数回の調整が発生することを抑制でき、制約係数の決定にかかる時間を短縮できる。
【0072】
【0073】
情報処理装置１００は、組合せ最適化問題を定式化することで、組合せ最適化問題に対応するエネルギー関数の情報を生成し、生成したエネルギー関数の情報を最適化装置２００に入力する。第２の実施の形態では、計算対象の組合せ最適化問題の例として、ＱＡＰが扱われる。
【0074】
情報処理装置１００は、ＣＰＵ１０１、ＲＡＭ１０２、ＨＤＤ１０３、ＩＯ（Input/Output）インタフェース１０４、ＧＰＵ１０５、入力インタフェース１０６、媒体リーダ１０７及びＮＩＣ（Network Interface Card）１０８を有する。ＣＰＵ１０１は、第１の実施の形態の処理部１２の一例である。ＲＡＭ１０２は、第１の実施の形態の記憶部１１の一例である。
【0075】
ＣＰＵ１０１は、プログラムの命令を実行するプロセッサである。ＣＰＵ１０１は、ＨＤＤ１０３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ１０２にロードし、プログラムを実行する。なお、ＣＰＵ１０１は複数のプロセッサコアを含んでもよい。また、情報処理装置１００は複数のプロセッサを有してもよい。複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」ということがある。
【0076】
ＲＡＭ１０２は、ＣＰＵ１０１が実行するプログラムやＣＰＵ１０１が演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。ＲＡＭ１０２は、たとえば、ＤＲＡＭである。なお、情報処理装置１００は、ＲＡＭ１０２以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。
【0077】
ＨＤＤ１０３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、及び、データを記憶する不揮発性の記憶装置である。なお、情報処理装置１００は、フラッシュメモリやＳＳＤ（Solid State Drive）などの他の種類の記憶装置を備えてもよく、複数の不揮発性の記憶装置を備えてもよい。
【0078】
ＩＯインタフェース１０４は、最適化装置２００と接続され、ＣＰＵ１０１からの命令に従って、最適化装置２００に対するデータの入出力を行う。たとえば、ＩＯインタフェース１０４は、ＣＰＵ１０１の命令に応じて、ＲＡＭ１０２のデータを最適化装置２００のレジスタまたはメモリに書き込んだり、最適化装置２００からデータを読み出して、ＲＡＭ１０２に書き込んだりする。ＩＯインタフェース１０４としては、たとえば、ＰＣＩ－ｅ（Peripheral Component Interconnect - Express）などが用いられる。
【0079】
ＧＰＵ１０５は、ＣＰＵ１０１からの命令に従って、情報処理装置１００に接続されたディスプレイ１１１に画像を出力する。ディスプレイ１１１としては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、プラズマディスプレイ、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなど、任意の種類のディスプレイを用いることができる。
【0080】
入力インタフェース１０６は、情報処理装置１００に接続された入力デバイス１１２から入力信号を取得し、ＣＰＵ１０１に出力する。入力デバイス１１２としては、マウス・タッチパネル・タッチパッド・トラックボールなどのポインティングデバイス、キーボード、リモートコントローラ、ボタンスイッチなどを用いることができる。また、情報処理装置１００に、複数の種類の入力デバイスが接続されていてもよい。
【0081】
媒体リーダ１０７は、記録媒体１１３に記録されたプログラムやデータを読み取る読み取り装置である。記録媒体１１３として、たとえば、磁気ディスク、光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。磁気ディスクには、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤが含まれる。光ディスクには、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）が含まれる。
【0082】
媒体リーダ１０７は、たとえば、記録媒体１１３から読み取ったプログラムやデータを、ＲＡＭ１０２やＨＤＤ１０３などの他の記録媒体にコピーする。読み取られたプログラムは、たとえば、ＣＰＵ１０１によって実行される。なお、記録媒体１１３は可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体１１３やＨＤＤ１０３を、コンピュータ読み取り可能な記録媒体ということがある。
【0083】
ＮＩＣ１０８は、ネットワーク３０に接続され、ネットワーク３０を介して他のコンピュータと通信を行うインタフェースである。ＮＩＣ１０８は、たとえば、スイッチやルータなどの通信装置とケーブルで接続される。
【0084】
最適化装置２００は、エネルギー関数の情報に基づいて、ＳＡ法やレプリカ交換法などのＭＣＭＣ法による基底状態探索をハードウェアにより行うアクセラレータである。最適化装置２００は、イジングマシンやボルツマンマシンなどと呼ばれてもよい。最適化装置２００は、ＦＰＧＡ２０１及びＲＡＭ２０２を有する。ＦＰＧＡ２０１は、イジングモデルで表される問題に対する解探索機能を実現する。最適化装置２００による解探索機能は、ＡＳＩＣやＧＰＵなどの他の種類の集積回路により実現されてもよい。ＲＡＭ２０２は、ＦＰＧＡ２０１の解探索に用いられる各種のデータを記憶するメモリである。ＲＡＭ２０２は、たとえば、ＳＲＡＭ（Static RAM）である。
【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】
ステップＳ５の処理では、制御部１３０は、調整係数の決定を行う。ステップＳ５の処理については後述する。ステップＳ５の処理後、制御部１３０は、ステップＳ６の処理を行う。
【0119】
ステップＳ６の処理では、制御部１３０は、Δ－Ｃ（ｘ）に調整係数を乗ずることで、制約係数αを算出する。
その後、制御部１３０は、式（１）をＱＵＢＯ形式で表したエネルギー関数の情報を生成する（ステップＳ７）。
その後、制御部１３０は、式（１）をＱＵＢＯ形式で表したエネルギー関数の情報を生成する（ステップＳ７）。
【0120】
制御部１３０は、ステップＳ７で生成したエネルギー関数の情報を探索部２１０に設定する（ステップＳ８）。
探索部２１０は、設定されたエネルギー関数の情報に基づいて、組合せ最適化演算を実行する。探索部２１０は、当該演算により、ＱＡＰの解を探索する（ステップＳ９）。
探索部２１０は、設定されたエネルギー関数の情報に基づいて、組合せ最適化演算を実行する。探索部２１０は、当該演算により、ＱＡＰの解を探索する（ステップＳ９）。
【0121】
探索部２１０は、探索により得られた解を制御部１３０に出力する（ステップＳ１０）。
制御部１３０は、探索部２１０から解を取得し、解出力部１４０を介して外部に出力し、処理を終了する。
制御部１３０は、探索部２１０から解を取得し、解出力部１４０を介して外部に出力し、処理を終了する。
【0122】
（比較例）
図９は、比較例の情報処理システムの処理の流れを示すフローチャートである。
図９には、従来のように、イジングマシンによる組合せ最適化問題の解の探索結果に基づいて、制約係数を調整する処理を繰り返すことで、制約係数を最適化していく手法を用いた処理の手順が示されている。
図９は、比較例の情報処理システムの処理の流れを示すフローチャートである。
図９には、従来のように、イジングマシンによる組合せ最適化問題の解の探索結果に基づいて、制約係数を調整する処理を繰り返すことで、制約係数を最適化していく手法を用いた処理の手順が示されている。
【0123】
ステップＳ２０，Ｓ２１の処理は、図６に示したステップＳ１，Ｓ２の処理と同じである。
ステップＳ２２の処理では、制約係数αの初期値（仮値）が設定される。その後、図６に示したステップＳ７～Ｓ９と同様の処理が行われ（ステップＳ２３，Ｓ２４，Ｓ２５）、組合せ最適化演算によって得られた割当行列ｘが出力される（ステップＳ２６）。
ステップＳ２２の処理では、制約係数αの初期値（仮値）が設定される。その後、図６に示したステップＳ７～Ｓ９と同様の処理が行われ（ステップＳ２３，Ｓ２４，Ｓ２５）、組合せ最適化演算によって得られた割当行列ｘが出力される（ステップＳ２６）。
【0124】
そして、割当行列ｘに基づいて、ＱＡＰの制約条件が満たされているか否かが判定される（ステップＳ２７）。制約条件が満たされていると判定された場合、ステップＳ２８の処理が行われ、制約条件が満たされていないと判定された場合、ステップＳ３０の処理が行われる。
【0125】
ステップＳ２８の処理では、組合せ最適化演算によって得られたコスト項の値が出力される。
その後、コスト項の値が所定の条件を満たしているか否かが判定され（ステップＳ２９）、条件が満たされていると判定された場合、ステップＳ３１の処理が行われ、条件が満たされていないと判定された場合、ステップＳ３０の処理が行われる。
その後、コスト項の値が所定の条件を満たしているか否かが判定され（ステップＳ２９）、条件が満たされていると判定された場合、ステップＳ３１の処理が行われ、条件が満たされていないと判定された場合、ステップＳ３０の処理が行われる。
【0126】
ステップＳ３０の処理では、制約係数αの更新が行われる。たとえば、制約充足解が得られていなければ制約係数を大きくするような調整が行われ、コスト項の値が所定の閾値よりも大きいには制約係数を小さくするような調整が行われる。ステップＳ３０の処理後、ステップＳ２３からの処理が繰り返される。
【0127】
ステップＳ３１の処理は、図６に示したステップＳ１０と同様の処理が行われる。
次に上記のような比較例の手法を用いて、ＱＡＰの各種のインスタンスを計算した例を示す。
次に上記のような比較例の手法を用いて、ＱＡＰの各種のインスタンスを計算した例を示す。
【0128】
図１０は、比較例における制約係数の値とインスタンスの計算結果との関係を示す図である。図１０には、レプリカ交換法を行うイジングマシンを用いて、問題サイズｎがｎ＝１２の５つのインスタンスに対して、制約係数の値を変えて計算を行った例が示されている。５つのインスタンスは、ＱＡＰＬＩＢから選択されたｔａｉ１２ａ、ｓｃｒ１２、ｒｏｕ１２、ｎｕｇ１２、ｈａｄ１２である。横軸は制約係数の値を表し、縦軸は１２８個のレプリカのうち制約充足解が得られたレプリカ数を表す。
【0129】
制約係数を小さい値から大きい値に変化させると、どのインスタンスに対する計算結果でも、制約係数が小さい部分では制約違反解のみしか得られていないが、制約係数を大きくしていくと制約充足解が得られるようになる。最終的には、ほぼ全てのレプリカで制約充足解が得られている。このように、制約係数を十分に大きくすれば、ほぼ確実に制約充足解を得ることができるが、一般的には、制約係数が小さいものほどコスト項の値は最適解に近い解（準最適解：コスト項の値が最適解の値に十分に近い解）が得られる。そのため、制約係数は、制約充足解が得られ始める値から、ほとんど制約充足解で一部制約違反解が残る程度の値（グラフの傾きが正の部分）が適切な値となる。図１０に示されているように、同じ問題サイズでも、インスタンスによって適切な制約係数の値は百未満から数万まで大きく異なることがわかる。
【0130】
たとえば、ｒｏｕ１２に対する求解の初期状態（図１０に示すような関係を計算する前の状態）では、制約係数の概数も不明であるため、十分大きそうな値を制約係数の初期値として設定することになる。制約係数の初期値を適当に、たとえば（１）５０，０００として組合せ最適化演算を行うと、１２８のレプリカの全てにおいて制約充足解が得られ、５０，０００は大きすぎるであろうことがわかる。次に、たとえば５０，０００を半分に減少して（２）２５，０００を制約係数の値として組合せ最適化演算を行う（図９のステップＳ２３の処理に戻る部分に対応している）と、１２８のレプリカ全てが制約違反解となる。このため、適切な制約係数（α）の値の範囲は、２５，０００＜α＜５０，０００であることがわかる。そこで次は、たとえば、２５，０００と５０，０００の間の（３）４０，０００を制約係数の値とし、組合せ最適化演算が行われる。しかしその場合、全てのレプリカで制約充足解が得られるため、やはり大きすぎることがわかる。続いて（４）３０，０００を制約係数の値として、組合せ最適化演算が行われる。この場合、一部のレプリカで制約充足解が得られ、制約係数は適切な値に近いことがわかる。この後、制約係数の値として、３０，０００と４０，０００の間や３０，０００より少し小さな値など何点かの値を用いた組合せ最適化演算が実行され、最終的に、たとえば、２９，５００が適切な制約係数として決定される。以上の操作はアルゴリズムとして組み込むことも可能である。
【0131】
一方、ｎｕｇ１２などを計算する場合、制約係数の初期値を５０，０００とすると、たとえば、５０，０００、２０，０００、１０，０００、５，０００、２，０００、１，０００、５００、２００、１００、５０まで変化させたときに制約違反解が得られる。このため、無駄な試行が非常に多くなる。
【0132】
人手ではなく、ベイズ推定などを用いて制約係数を最適化することも考えられるが、多数回の試行を要する点は変わらない。
上記のような比較例に対して、図６に示したような処理を行う本実施の形態の情報処理システム２０では、計算したΔ－Ｃ（ｘ）に基づいて制約係数αを決定することで、制約係数αの調整時間を短縮できる。なお、図６には、調整係数を決定するステップＳ５が含まれているが、以下に説明するように、適切な調整係数の探索範囲は限定されるため、ステップＳ５の処理による計算時間の増加は、比較例の場合と比べて少ない。
上記のような比較例に対して、図６に示したような処理を行う本実施の形態の情報処理システム２０では、計算したΔ－Ｃ（ｘ）に基づいて制約係数αを決定することで、制約係数αの調整時間を短縮できる。なお、図６には、調整係数を決定するステップＳ５が含まれているが、以下に説明するように、適切な調整係数の探索範囲は限定されるため、ステップＳ５の処理による計算時間の増加は、比較例の場合と比べて少ない。
【0133】
（調整係数の決定）
次に、調整係数の決定について説明する。
イジングマシンで求解が行われるインスタンス（ＱＵＢＯ行列を用いて表される）は、課題によって「統計的性質」(行列要素の平均、分散、行列の疎密、値の分布など)が異なる。このため、Δ－Ｃ（ｘ）をインスタンス系列に合わせて調整することで制約係数αを決定することが望ましい。
次に、調整係数の決定について説明する。
イジングマシンで求解が行われるインスタンス（ＱＵＢＯ行列を用いて表される）は、課題によって「統計的性質」(行列要素の平均、分散、行列の疎密、値の分布など)が異なる。このため、Δ－Ｃ（ｘ）をインスタンス系列に合わせて調整することで制約係数αを決定することが望ましい。
【0134】
ＱＡＰとして定式化される課題には、工場用地への工場の割り当てのみならず、様々なものがある。たとえば、病院における病室と施設の配置や、空港での乗り換え乗客の延べ移動距離を最小化するゲート配置、プリント板の配置配線問題、キーボードのキー配列、テーブルへの着座配置、タービンブレードの重量バランスなどがある。これらの課題のそれぞれにおいて「統計的性質」が異なる。
【0135】
なお、ＱＡＰＬＩＢはそのような課題を抽象化し、課題ごとにサイズの異なるシリーズとして作成されたライブラリで、解精度や求解速度のベンチマークに用いられている。
図１１は、本実施の形態における調整係数の値とインスタンスの計算結果との関係を示す図である。図１１には、第２の実施の形態の情報処理システム２０により、問題サイズｎがｎ＝１２の５つのインスタンスに対して、Ｃ（ｘ）の減少量の期待値Δ－Ｃ（ｘ）を求めた後、調整係数の値を変えてレプリカ交換法による計算を行った例が示されている。５つのインスタンスは、ＱＡＰＬＩＢから選択されたｔａｉ１２ａ、ｓｃｒ１２、ｒｏｕ１２、ｎｕｇ１２、ｈａｄ１２である。横軸は調整係数の値を表し、縦軸は１２８個のレプリカのうち制約充足解が得られたレプリカ数を表す。
図１１は、本実施の形態における調整係数の値とインスタンスの計算結果との関係を示す図である。図１１には、第２の実施の形態の情報処理システム２０により、問題サイズｎがｎ＝１２の５つのインスタンスに対して、Ｃ（ｘ）の減少量の期待値Δ－Ｃ（ｘ）を求めた後、調整係数の値を変えてレプリカ交換法による計算を行った例が示されている。５つのインスタンスは、ＱＡＰＬＩＢから選択されたｔａｉ１２ａ、ｓｃｒ１２、ｒｏｕ１２、ｎｕｇ１２、ｈａｄ１２である。横軸は調整係数の値を表し、縦軸は１２８個のレプリカのうち制約充足解が得られたレプリカ数を表す。
【0136】
図１０に示した比較例の場合と比べて、インスタンス間での適切な調整係数の値のばらつきは大幅に小さい。ただ、詳細に見ると、それぞれのインスタンスで統計的性質が異なるため、適切な調整係数（グラフの傾きが正の範囲）には数十パーセント程度の違いがある。このため、第２の実施の形態の情報処理システム２０がより良好な解を得るためには、インスタンスごとに適切な調整係数を定めることが好ましい。
【0137】
図１２は、調整係数の決定処理の一例の流れを示すフローチャートである。
まず、調整係数の初期値、増加量及び上限の設定が行われる（ステップＳ４０）。ステップＳ４０の処理では、たとえば、ユーザにより入力された初期値、増加量または上限に基づいて設定が行われてもよいし、予め記憶部１２０に記憶された初期値、増加量または上限に基づいて設定が行われてもよい。
まず、調整係数の初期値、増加量及び上限の設定が行われる（ステップＳ４０）。ステップＳ４０の処理では、たとえば、ユーザにより入力された初期値、増加量または上限に基づいて設定が行われてもよいし、予め記憶部１２０に記憶された初期値、増加量または上限に基づいて設定が行われてもよい。
【0138】
その後、制御部１３０と探索部２１０により、図６に示したステップＳ６～Ｓ９の処理と同様のステップＳ４１～Ｓ４４の処理が行われ、イジングモデルの基底状態の探索結果である割当行列ｘが探索部２１０から出力される（ステップＳ４５）。
【0139】
ステップＳ４５の処理の後、制御部１３０は、調整係数が上限に達したか否かを判定する（ステップＳ４６）。
制御部１３０は、調整係数が上限に達していないと判定した場合には、調整係数を増加量分、増加させ（ステップＳ４７）、ステップＳ４１からの処理を繰り返す。
制御部１３０は、調整係数が上限に達していないと判定した場合には、調整係数を増加量分、増加させ（ステップＳ４７）、ステップＳ４１からの処理を繰り返す。
【0140】
制御部１３０は、調整係数が上限に達したと判定した場合には、これまでに得られた探索結果である割当行列ｘに基づいて、最適な調整係数を決定し（ステップＳ４８）、調整係数の決定処理を終了する。制御部１３０は、これまでに得られた割当行列ｘに基づいて、たとえば、図１１に示したように、制約充足解が得られたレプリカ数が調整係数の変化に応じて増加する範囲（傾きが正となる範囲）で調整係数を決定することで、適切な調整係数が得られる。
【0141】
適切な調整係数の探索範囲は、１×１００前後に限られ、しかも探索の刻み（上記の調整係数の増加量）は０．０５～０．１程度で十分である。また、調整係数の調整回数は、ステップＳ４０の処理で設定された初期値、増加量及び上限によって決まっている。
【0142】
このため、比較例において制約係数αを更新する回数よりも、調整係数の調整回数を大幅に少なくできることが期待できる。
さらに、一度、調整係数を定めれば同種のインスタンスの組合せ最適化演算の実行時には、その調整係数を用いることができるため、調整係数決定のための試行は初回のみ行えばよい。
さらに、一度、調整係数を定めれば同種のインスタンスの組合せ最適化演算の実行時には、その調整係数を用いることができるため、調整係数決定のための試行は初回のみ行えばよい。
【0143】
図１３は、調整係数の値と問題サイズの異なる同種のインスタンスの計算結果との関係を示す図である。図１３には、第２の実施の形態の情報処理システム２０により、問題サイズの異なる１２のｔａｉＮａ（Ｎは問題サイズ）に対して、Ｃ（ｘ）の減少量の期待値Δ－Ｃ（ｘ）を求めた後、調整係数の値を変えてレプリカ交換法による計算を行った例が示されている。横軸は調整係数の値を表し、縦軸は１２８個のレプリカのうち制約充足解が得られたレプリカ数を表す。なお、図１３には、各インスタンスについての計算結果が、縦軸方向に１５０ずつずらして表示されている。
【0144】
図１３の例では、たとえば、問題サイズの一番小さいｔａｉ１０ａにおいて、制約充足解が得られ始める程度の調整係数（１．３程度）を用いることで、１２の異なる問題サイズの同種のインスタンス（ｔａｉ）に対して、制約充足解を得ることができる。
【0145】
このように、調整係数を求める際は計算対象の問題の問題サイズを縮小して適切な調整係数を探索し、その調整係数を用いて実際の問題の求解を行うことにより、調整係数探索時間を大幅に削減することが可能となる。
【0146】
調整係数を求める際には、最適化しようとする問題を縮小したモデル（たとえば工場の用地への割り当てでは工場数、用地数を減らしたもの）を用いることが好ましい。ただ、抽象的な問題でもフロー行列、距離行列に統計的に大きな偏りがなければ、それぞれの行列から適当な行または列を抜き出して構成した小サイズのフロー行列、距離行列を用いて調整係数を求めることも可能である。
【0147】
なお、上記の調整係数の決定例では、調整係数の増加量と上限を設定し、調整係数が上限に達するまで、調整係数を増加させる方法が用いられたが、調整係数の減少量と下限を設定し、調整係数が下限に達するまで、調整係数を減少させる方法を用いてもよい。
【0148】
（ＱＵＢＯ生成例）
次に、決定した制約係数αに基づいたＱＵＢＯの生成例を説明する。
制約項は、たとえば、前述の式（８）のように表せる。式（８）の右辺の１項目のαに乗ずる２次多項式は、以下の式（１４）のように表すことができる。
次に、決定した制約係数αに基づいたＱＵＢＯの生成例を説明する。
制約項は、たとえば、前述の式（８）のように表せる。式（８）の右辺の１項目のαに乗ずる２次多項式は、以下の式（１４）のように表すことができる。
【0149】
【数14】
【0150】
式（１４）において、右辺の２次項は以下の式（１５）のように表すことができる。
【0151】
【数15】
【0152】
式（１５）において２次項は、状態変数ベクトルの外積の上三角（ｉ＞ｊ）の非対角要素の範囲（後述の図１４参照）で状態変数の和をとり、下三角の分、２倍したものである。また、式（１５）において、ｘｉｋ
２は対角要素（後述の図１４参照）である。ｘｉｋは０または１のバイナリ変数であるため、ｘｉｋ
２＝ｘｉｋである。以上のことから、式（１４）は、以下の式（１６）のように変形できる。
【0153】
【数16】
【0154】
式（８）の右辺の２項目のβに乗ずる２次多項式も上記と同様の理由により、以下の式（１７）のように表せる。
【0155】
【数17】
【0156】
式（１６）における２次項に足し合わされるｘｉｋｘｊｋと、１次項とを行列で表すと以下のようになる。
図１４は、式（１６）における２次項に足し合わされるｘｉｋｘｊｋと、１次項とを行列で表した図である。図１４では問題サイズｎがｎ＝４である場合の例が示されている。
図１４は、式（１６）における２次項に足し合わされるｘｉｋｘｊｋと、１次項とを行列で表した図である。図１４では問題サイズｎがｎ＝４である場合の例が示されている。
【0157】
領域４０内の２状態変数の積が２次項に足し合わされる。図１４では、２次項に足し合わされる２状態変数の積を表す行列要素が、ｋ（＝１～４）の値ごとに異なるハッチングで表されている。
【0158】
図１４では、１次項は、行列の対角要素として、ｋ（＝１～４）の値ごとに異なるハッチングで表されている。
式（１７）における２次項に足し合わされるｘｉｋｘｉｌと、１次項とを行列で表すと以下のようになる。
式（１７）における２次項に足し合わされるｘｉｋｘｉｌと、１次項とを行列で表すと以下のようになる。
【0159】
図１５は、式（１７）における２次項に足し合わされるｘｉｋｘｉｌと、１次項とを行列で表した図である。図１５でも問題サイズｎがｎ＝４である場合の例が示されている。
領域４１内の２状態変数の積が２次項に足し合わされる。図１５では、２次項に足し合わされる２状態変数の積を表す行列要素が、ｉ（＝１～４）の値ごとに異なるハッチングで表されている。
領域４１内の２状態変数の積が２次項に足し合わされる。図１５では、２次項に足し合わされる２状態変数の積を表す行列要素が、ｉ（＝１～４）の値ごとに異なるハッチングで表されている。
【0160】
【0161】
図１４、図１５の結果を組み合わせると、制約項のＱＵＢＯ行列は、図１６のように表される。制約項についてのＱＵＢＯ行列は、対角要素は－（α＋β）であり、上三角部分（ｉ＞ｊ）の各小行列の対角要素は、２αである。また、対角の４つの小行列のそれぞれにおける上三角部分の各要素は２βである。なお、前述のようにα＝βと考えてよい。
【0162】
図１７は、エネルギー関数全体のＱＵＢＯ行列の例を示す図である。
図１７のエネルギー関数全体のＱＵＢＯ行列は、図７に示したコスト項のＱＵＢＯ行列と、図１６に示した制約項のＱＵＢＯ行列とを足し合わせることで得られる。
図１７のエネルギー関数全体のＱＵＢＯ行列は、図７に示したコスト項のＱＵＢＯ行列と、図１６に示した制約項のＱＵＢＯ行列とを足し合わせることで得られる。
【0163】
制御部１３０は、このようなＱＵＢＯ行列をエネルギー関数の情報として出力し、探索部２１０に設定する。
以上のような第２の実施の形態の情報処理システム２０によれば、第１の実施の形態の情報処理システム１０と同様の効果が得られる。すなわち、適切な制約係数を得るために多数回の調整が発生することを抑制でき、制約係数の決定にかかる時間を短縮できる。
以上のような第２の実施の形態の情報処理システム２０によれば、第１の実施の形態の情報処理システム１０と同様の効果が得られる。すなわち、適切な制約係数を得るために多数回の調整が発生することを抑制でき、制約係数の決定にかかる時間を短縮できる。
【0164】
（第３の実施の形態）
図１８は、第３の実施の形態の情報処理システムの一例の処理の流れを示すフローチャートである。
図１８は、第３の実施の形態の情報処理システムの一例の処理の流れを示すフローチャートである。
【0165】
ステップＳ５０～Ｓ５１の処理は、図６に示したステップＳ１，Ｓ２の処理と同じである。
ステップＳ５１の処理後、制御部１３０は、行列変調を行う（ステップＳ５２）。ステップＳ５２の処理では、フロー行列Ｆと距離行列Ｄのクロネッカー積によって生成されたＱｉｊの変調が行われる。変調の例については後述する。
ステップＳ５１の処理後、制御部１３０は、行列変調を行う（ステップＳ５２）。ステップＳ５２の処理では、フロー行列Ｆと距離行列Ｄのクロネッカー積によって生成されたＱｉｊの変調が行われる。変調の例については後述する。
【0166】
ステップＳ５２の処理後、制御部１３０は、変調後のＱｉｊに基づいて、図６に示したステップＳ３の処理と同様に、Δ－Ｃ（ｘ）を計算する（ステップＳ５３）。第３の実施の形態の情報処理システムでは、後述の理由によって調整係数を１．０とすることができ、Δ－Ｃ（ｘ）を制約係数αとすることができる。
【0167】
【0168】
【数18】
【0169】
組合せ最適化問題の計算では、常に同じ定数Ｃをコスト項から加減しても最適解（ＱＡＰの場合は最適な割り当て）は変わらず、最適解におけるコスト項の値は、コスト項に加減した定数Ｃの分だけの差が生じるだけである。さらに、最適解に限らず、全ての解とそれらの順位（大小関係）は当然に保たれる。
【0170】
ところで、ＱＡＰでは、前述の式（９）のように、フロー行列Ｆと距離行列Ｄのクロネッカー積によって表される行列Ｑは、ｎ×ｎ個の小行列ｆｉｊＤからなる。ｉ，ｊで識別される任意の２つの割当要素（工場など）を考えたとき、割り当て可能な割当先（用地など）は、２つの割当先の全ての組合せである。小行列ｆｉｊＤは、ｉ，ｊで識別される任意の２つの割当要素に対する２つの割当先の全ての組合せに対応したコスト項を表す。
【0171】
ＱＡＰの制約充足解では、それぞれの小行列ｆｉｊＤから対角要素以外の何れか１つの要素ｆｉｊｄｋｌが選ばれる。すなわち制約充足解のコスト項には、それぞれの小行列ｆｉｊＤから対角要素以外の何れか１つの要素ｆｉｊｄｋｌの値のみが加算される。
【0172】
このため、制約充足解については、その解を構成する全てのｆｉｊｄｋｌから同じ定数を加減しても、上記のように単純にコスト項から常に定数Ｃを加減した場合と同様に、同一の最適解となる。また、全ての解とそれらの順位も保たれる。
【0173】
したがって、行列Ｑにおいて、前述の有効な要素からなるｎ（ｎ－１）／２個の小行列ｆｉｊＤごとに、それぞれの要素ｆｉｊｄｋｌから任意の同一定数を加減してもしなくても、制約充足解については同一の解（割り当て）及びコスト項順位となる。このような、それぞれの要素ｆｉｊｄｋｌから任意の同一定数を加減する操作が、行列Ｑの変調に相当する。
【0174】
一方、フロー行列Ｆと距離行列Ｄの順番を入れ替えてクロネッカー積を実施した場合に得られる小行列ｄｋｌＦは、ｋ，ｌで識別される任意の２つの割当先に対する２つの割当要素の全ての組合せに対応したコスト項を表す。しかし、クロネッカー積は可換ではないため、小行列ｄｋｌＦの要素は、式（９）に示した行列Ｑにおいては小行列を構成せず、行列Ｑの各小行列ｆｉｊＤに、１要素ずつ配分されている。
【0175】
【0176】
たとえば、ｄ２ｎＦの要素ｆ１２ｄ２ｎは小行列ｆ１２Ｄに配分されており、ｄ２ｎＦの要素ｆ１ｎｄ２ｎは小行列ｆ１ｎＤに配分されている。
ＱＡＰの制約充足解では、小行列ｆｉｊＤの場合と同様に、小行列ｄｋｌＦから対角要素以外の何れか１つの要素ｆｉｊｄｋｌが選ばれる。すなわち制約充足解のコスト項には、小行列ｄｋｌＦの対角要素以外の何れか１つの要素ｆｉｊｄｋｌの値が加算される。このため、制約充足解については、その解を構成する全てのｆｉｊｄｋｌから同じ定数を加減しても、上記のように単純にコスト項から常に定数Ｃを加減した場合と同様に、同一の最適解となる。また、全ての解とそれらの順位も保たれる。なお、小行列ｆｉｊＤの下三角の要素は、上三角の要素と同じであるため用いなくてもよい（無効な要素とすることができる）。
ＱＡＰの制約充足解では、小行列ｆｉｊＤの場合と同様に、小行列ｄｋｌＦから対角要素以外の何れか１つの要素ｆｉｊｄｋｌが選ばれる。すなわち制約充足解のコスト項には、小行列ｄｋｌＦの対角要素以外の何れか１つの要素ｆｉｊｄｋｌの値が加算される。このため、制約充足解については、その解を構成する全てのｆｉｊｄｋｌから同じ定数を加減しても、上記のように単純にコスト項から常に定数Ｃを加減した場合と同様に、同一の最適解となる。また、全ての解とそれらの順位も保たれる。なお、小行列ｆｉｊＤの下三角の要素は、上三角の要素と同じであるため用いなくてもよい（無効な要素とすることができる）。
【0177】
したがって、小行列ｆｉｊＤと同様に、有効な要素からなるｎ（ｎ－１）／２個の小行列ｄｋｌＦごとに、それぞれの要素ｆｉｊｄｋｌから任意の同一定数を加減してもしなくても、制約充足解については同一の解（割り当て）及びコスト項順位となる。このような、ｄｋｌＦのそれぞれの要素ｆｉｊｄｋｌから任意の同一定数を加減する操作も、行列Ｑの変調に相当する。
【0178】
これらの行列変調によりｆｉｊｄｋｌの値の偏差を縮減して解空間を滑らかにすることで、制約項の値を削減して解の探索効率を向上させることができる。これは、統計的な観点から以下のようにいうことができる。元々インスタンス系列ごとに異なる偏差を持っていたｆｉｊｄｋｌが、行列変調で偏差の縮減を行うことで、元のインスタンス系列によらず互いに統計的性質が類似し、かつ元のインスタンスと同一の解をもつ別のインスタンス系列に変換されている。
【0179】
このように、制約係数の決定に際して、計算対象のインスタンスに行列変調を施すことにより統計的性質を揃えることができる。これにより、第２の実施の形態に用いていた、インスタンス系列による統計的な違いを補正するために用いた調整係数を、インスタンスの系列依存性を軽減しつつ、ほとんど一定の値とすることが可能である。
【0180】
図２０は、行列変調の一例の流れを示すフローチャートである。
制御部１３０は、フロー行列Ｆと距離行列Ｄのクロネッカー積により得られる行列Ｑの要素の平均値を算出する（ステップＳ６０）。また、制御部１３０は、各小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値を算出する（ステップＳ６１）。そして制御部１３０は、各小行列ｆｉｊＤ（ｉ＜ｊ）について、（小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値－行列Ｑの要素の平均値）を、小行列ｆｉｊＤ（ｉ＜ｊ）の各要素から差し引く（ステップＳ６２）。なお、制御部１３０は、（小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値－行列Ｑの要素の平均値）が負の場合には、この差分値を、小行列ｆｉｊＤ（ｉ＜ｊ）の各要素に加える。
制御部１３０は、フロー行列Ｆと距離行列Ｄのクロネッカー積により得られる行列Ｑの要素の平均値を算出する（ステップＳ６０）。また、制御部１３０は、各小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値を算出する（ステップＳ６１）。そして制御部１３０は、各小行列ｆｉｊＤ（ｉ＜ｊ）について、（小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値－行列Ｑの要素の平均値）を、小行列ｆｉｊＤ（ｉ＜ｊ）の各要素から差し引く（ステップＳ６２）。なお、制御部１３０は、（小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値－行列Ｑの要素の平均値）が負の場合には、この差分値を、小行列ｆｉｊＤ（ｉ＜ｊ）の各要素に加える。
【0181】
次に制御部１３０は、各小行列ｄｋｌＦ（ｋ＜ｌ）の要素の平均値を算出する（ステップＳ６３）。そして制御部１３０は、各小行列ｄｋｌＦ（ｋ＜ｌ）について、（小行列ｄｋｌＦ（ｋ＜ｌ）の要素の平均値－行列Ｑの要素の平均値）を、小行列ｄｋｌＦ（ｋ＜ｌ）の各要素から差し引く（ステップＳ６４）。なお、制御部１３０は、（小行列ｄｋｌＦ（ｋ＜ｌ）の要素の平均値－行列Ｑの要素の平均値）が負の場合には、この差分値を、小行列ｄｋｌＦ（ｋ＜ｌ）の各要素に加える。
【0182】
最後に、制御部１３０は、行列Ｑにおける最小値をＱの全要素から差し引く（ステップＳ６５）。
上記の行列変調操作において、ステップＳ６２の処理は、各小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値が行列Ｑ全体の要素の平均値に一致するような値を各小行列の各要素から加減することに相当する。また、ステップＳ６４の処理は、各小行列ｄｋｌＦ（ｋ＜ｌ）の要素の平均値が行列Ｑ全体の要素の平均値に一致するような値を各小行列ｄｋｌＦ（ｋ＜ｌ）から加算または減算することに相当する。
上記の行列変調操作において、ステップＳ６２の処理は、各小行列ｆｉｊＤ（ｉ＜ｊ）の要素の平均値が行列Ｑ全体の要素の平均値に一致するような値を各小行列の各要素から加減することに相当する。また、ステップＳ６４の処理は、各小行列ｄｋｌＦ（ｋ＜ｌ）の要素の平均値が行列Ｑ全体の要素の平均値に一致するような値を各小行列ｄｋｌＦ（ｋ＜ｌ）から加算または減算することに相当する。
【0183】
なお、上記の処理の順序は一例であり、適宜処理の順序を変えてもよい。
図２１は、行列変調を用いた場合の調整係数の値とインスタンスの計算結果との関係を示す図である。図２１には、上記の方法により行列変調を行った行列Ｑに基づいてＣ（ｘ）の減少量の期待値Δ－Ｃ（ｘ）を求めた後、調整係数の値を変えてレプリカ交換法による計算を行った例が示されている。計算対象のインスタンスは、ＱＡＰＬＩＢから選択されたｔａｉ１２ａ、ｓｃｒ１２、ｒｏｕ１２、ｎｕｇ１２、ｈａｄ１２である。横軸は調整係数の値を表し、縦軸は１２８個のレプリカのうち制約充足解が得られたレプリカ数を表す。
図２１は、行列変調を用いた場合の調整係数の値とインスタンスの計算結果との関係を示す図である。図２１には、上記の方法により行列変調を行った行列Ｑに基づいてＣ（ｘ）の減少量の期待値Δ－Ｃ（ｘ）を求めた後、調整係数の値を変えてレプリカ交換法による計算を行った例が示されている。計算対象のインスタンスは、ＱＡＰＬＩＢから選択されたｔａｉ１２ａ、ｓｃｒ１２、ｒｏｕ１２、ｎｕｇ１２、ｈａｄ１２である。横軸は調整係数の値を表し、縦軸は１２８個のレプリカのうち制約充足解が得られたレプリカ数を表す。
【0184】
上記の方法により行列変調を行った行列Ｑを用いた場合、図２１に示されているように、５つのインスタンスについて調整係数をほぼ１．０にすることが可能であった。つまり、図１８に示したように、制約係数αの決定のための試行錯誤をしなくても、Δ－Ｃ（ｘ）をそのまま制約係数αとすることができる。
【0185】
ただし、調整係数が１．０になるのは図２０に示したような行列変調操作が行われる場合であって、他の行列変調操作が行われる場合は、調整係数が１．０とは異なる値となる場合もありうる。
【0186】
また、制約係数αを無調整で解を得る前に、あるいは解を得た後に再度、制約係数αを調整してさらによいコスト項の値を目指すことも、もちろん可能である。
なお、第１の実施の形態の情報処理は、処理部１２にプログラムを実行させることで実現できる。また、第２及び第３の実施の形態の情報処理は、ＣＰＵ１０１にプログラムを実行させることで実現できる。プログラムは、コンピュータ読み取り可能な記録媒体１１３に記録できる。
なお、第１の実施の形態の情報処理は、処理部１２にプログラムを実行させることで実現できる。また、第２及び第３の実施の形態の情報処理は、ＣＰＵ１０１にプログラムを実行させることで実現できる。プログラムは、コンピュータ読み取り可能な記録媒体１１３に記録できる。
【0187】
たとえば、プログラムを記録した記録媒体１１３を配布することで、プログラムを流通させることができる。また、プログラムを他のコンピュータに格納しておき、ネットワーク経由でプログラムを配布してもよい。コンピュータは、たとえば、記録媒体１１３に記録されたプログラムまたは他のコンピュータから受信したプログラムを、ＲＡＭ１０２やＨＤＤ１０３などの記憶装置に格納し（インストールし）、当該記憶装置からプログラムを読み込んで実行してもよい。
【0188】
以上、実施の形態に基づき、本発明の情報処理システム、情報処理方法及びプログラムの一観点について説明してきたが、これらは一例にすぎず、上記の記載に限定されるものではない。
【符号の説明】
【0189】
１０ 情報処理システム
１１ 記憶部
１２ 処理部
１３ 探索部
１１ 記憶部
１２ 処理部
１３ 探索部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】