特許 7536710 求解装置、求解方法およびプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-08-09

(45)【発行日】2024-08-20

(54)【発明の名称】求解装置、求解方法およびプログラム

(51)【国際特許分類】

   G06N 99/00 20190101AFI20240813BHJP

   G06F 17/10 20060101ALI20240813BHJP

【ＦＩ】

G06N99/00 180

G06F17/10 Z

【請求項の数】 16

(21)【出願番号】P 2021087125

(22)【出願日】2021-05-24

(65)【公開番号】P2022180174

(43)【公開日】2022-12-06

【審査請求日】2023-11-15

(73)【特許権者】

【識別番号】000003078

【氏名又は名称】株式会社東芝

(73)【特許権者】

【識別番号】301063496

【氏名又は名称】東芝デジタルソリューションズ株式会社

(74)【代理人】

【識別番号】110002147

【氏名又は名称】弁理士法人酒井国際特許事務所

(72)【発明者】

【氏名】鈴木 賢

【審査官】福西 章人

(56)【参考文献】

【文献】国際公開第２０２０／１９６８６２（ＷＯ，Ａ１）

【文献】後藤 隼人，シミュレーテッド分岐アルゴリズム，日本オペレーションズ・リサーチ学会 ２０２０年春季研究発表会アブストラクト集，2020年，pp.298-299，ISSN 1883-1893

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０６Ｎ ３／００－９９／００

Ｇ０６Ｆ １７／１０

(57)【特許請求の範囲】

【請求項1】

０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く求解装置であって、
第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する出力部と、
を備え、
前記複数の要素は、前記複数の離散変数に対応し、
前記第１変数および前記第２変数のそれぞれは、実数により表され、
前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、
前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新し、
前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
求解装置。

【請求項2】

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

【請求項3】

前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第１変数を更新した後に前記第２変数を更新し、
前記第２変数を更新した後に前記第２変数から前記補正値を減算する
請求項２に記載の求解装置。

【請求項4】

前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、
前記第２変数を更新した後に前記第１変数を更新し、
前記第２変数を更新した後、且つ、前記第１変数を更新する前に、前記第２変数から前記補正値を減算する
請求項２に記載の求解装置。

【請求項5】

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

【請求項6】

前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、
前記複数の要素のそれぞれの前記第１変数と、対象要素と前記複数の要素のそれぞれとの組毎に前記０－１組合せ最適化問題により定められる作用係数とに基づき、外力を算出し、
時間経過に従って増加する関数に基づき定まる値と、前記対象要素の前記第１変数とを乗算した時間発展値を算出し、
前記直前時刻における前記第２変数に、前記時間発展値と前記外力とを加算した値に前記単位時間を乗算した値を、加算することにより、前記対象時刻における前記第２変数を算出する
請求項５に記載の求解装置。

【請求項7】

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

【数1】

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

【請求項8】

前記０－１組合せ最適化問題は、０－１二次計画問題であり、
前記更新部は、前記ｉ番目の要素の前記対象時刻における前記第２変数を式（１０３）または式（１０４）により算出し、

【数2】

ｇ_ｉ（ｔ）およびｇ_ｉ（ｔ＋Δｔ）は、前記時間発展値であり、
ｆ_ｉ（ｔ＋Δｔ）は、式（１０５）により表され、ｆ_ｉ（ｔ）は、式（１０６）により表され、

【数3】

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

【数4】

ｊは、１からＮまでの整数であり、
ｈ_ｉは、Ｎ個の局所磁場係数を含む予め定められた配列に含まれる、ｉ番目の局所磁場係数であり、
Ｊ_ｉ，ｊは、Ｎ×Ｎ個の結合係数を含む予め定められた行列に含まれる、ｉ行、ｊ列の結合係数であり、
ｃは、予め定められた実数であり、
ｘ_ｊ（ｔ）は、前記Ｎ個の離散変数のうちのｊ番目の離散変数に対応するｊ番目の要素の前記直前時刻における前記第１変数であり、
ｘ_ｊ（ｔ＋Δｔ）は、前記ｊ番目の要素の前記対象時刻における前記第１変数であり、
α（ｔ）は、ｔを変数とする予め定められた関数である
請求項７に記載の求解装置。

【請求項9】

ｇ_ｉ（ｔ）は、式（１０９）により表され、
ｇ_ｉ（ｔ＋Δｔ）は、式（１１０）により表され、

【数5】

ｐ（ｔ）は、ｔを変数とする予め定められた関数であり、ｔに従って増加し、ｔが前記初期時刻において０となり、ｔが前記終了時刻において１となり、
Ｋは、予め定められた定数である
請求項８に記載の求解装置。

【請求項10】

前記１または複数の制約条件のうちのｍ番目の制約条件は、式（１１１）により表され、

【数6】

ｍは、１以上の整数であり、
ｘ_ｉは、ｉ番目の離散変数であり、
ω_ｍ，ｉは、ｉ番目の離散変数に乗算される係数であり、
Ｗ_ｍは、予め定められた定数である
請求項８または９に記載の求解装置。

【請求項11】

前記更新部は、前記ｉ番目の要素の前記第２変数に対する前記補正値として、式（１１２）により表されるｒ_ｍ，ｉを算出し、

【数7】

Ａ_ｍは、前記ｍ番目の制約条件に対して予め定められた係数であり、
ｋは、予め定められた１以上の整数であり、
Ｓ_ｍは、式（１１３）または式（１１４）により表される

【数8】

請求項１０に記載の求解装置。

【請求項12】

前記更新部は、前記ｉ番目の要素の前記第２変数に対する前記補正値として、式（１１５）により表されるｒ_ｍ，ｉを算出し、

【数9】

Ａ_ｍは、前記ｍ番目の制約条件に対して予め定められた係数であり、
ｋは、予め定められた１以上の整数であり、
Ｓ_ｍは、式（１１６）または式（１１７）により表され、

【数10】

Ｂ_ｍは、０以上１以下の値であり、前記ｍ番目の制約条件に対して予め定められた係数であり、
Ｅ_ｍは、式（１１８）により表される

【数11】

請求項１０に記載の求解装置。

【請求項13】

前記１または複数の制約条件のうちの一つの制約条件は、式（１１９）により表され、

【数12】

ｘ_ｉは、ｉ番目の離散変数であり、
ｘ_ｊは、ｊ番目の離散変数であり、
Ｑ_ｉ，ｊは、半正定値行列に含まれるｉ行、ｊ列の値であり、
ｑは、予め定められた定数である
請求項８または９に記載の求解装置。

【請求項14】

前記更新部は、前記ｉ番目の要素の前記第２変数に対する前記補正値として、式（１２０）または（１２１）により表されるｒ_ｉを算出し、

【数13】

Ａは、対応する制約条件に対して予め定められた係数であり、
ｋは、予め定められた１以上の整数であり、
Ｓは、式（１２０）の場合は式（１２２）により表され、式（１２１）の場合には式（１２３）により表される

【数14】

請求項１３に記載の求解装置。

【請求項15】

情報処理装置により、０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く求解方法であって、
前記情報処理装置により、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新ステップと、
前記情報処理装置により、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する出力ステップと、
を含み、
前記複数の要素は、前記複数の離散変数に対応し、
前記第１変数および前記第２変数のそれぞれは、実数により表され、
前記単位時間毎の更新ステップにおいて、前記情報処理装置は、
前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、
前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新し、
前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
求解方法。

【請求項16】

情報処理装置を、０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く求解装置として機能させるためのプログラムであって、
前記情報処理装置を、
第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する更新部と、
前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する出力部と、
して機能させ、
前記複数の要素は、前記複数の離散変数に対応し、
前記第１変数および前記第２変数のそれぞれは、実数により表され、
前記単位時間毎の更新処理において、前記更新部は、
前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、
前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新し、
前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する
プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

従来、一次不等式による制約の下でイジング問題を解く方法として、スラック変数を利用した定式化をする方法が知られている。この方法は、一次不等式の定数項をＷ_ｍとした場合、目的関数に、Ｗ個（Ｗは自然数）のスラック変数により表されたペナルティ項を加算して、イジング問題を解く。従って、この方法は、目的関数にＮ個の決定変数が含まれるならば、Ｎ＋Ｗ個のイジングスピンを含むイジング問題を解く。また、例えば、Ｍ個の一次不等式による制限の下でイジング問題を解くのであれば、この方法は、Ｎ＋｛Ｗ_１＋…＋Ｍ_ｍ＋…＋Ｗ_Ｍ｝個のイジングスピンを含むイジング問題を解く。なお、Ｍ_ｍは、Ｍ個の一次不等式のうちのｍ番目の一次不等式の定数項であり、自然数を表す。

【0003】

株式のポートフォリオ最適化問題は、リスク回避のために、銘柄の投資割合を決定変数とし、複数の決定変数のそれぞれに重みを乗算した合計を上限値以下とする制約条件の下で、解かれる場合が多い。制約条件は、下記の式のような一次不等式により表される。
ｘ_Ａ１・ω_Ａ１＋ｘ_Ａ２・ω_Ａ２＋…＋ｘ_Ｓｍ・ω_Ｓｍ＋…≦Ｗ_ｍ

【0004】

この式において、Ｗ_ｍは、上限値を表す。ｘ_Ａ１は、銘柄Ａの投資割合を例えば１％とする決定変数を表す。ｘ_Ａ２は、銘柄Ａの投資割合を例えば２％とする決定変数を表す。ｘ_Ｓｓは、銘柄Ｓの投資割合をｓ％とする決定変数を表す。また、ω_Ａ１、ω_Ａ２およびω_Ｓｍは、重みを表す。

【0005】

例えば、株式のポートフォリオ最適化問題は、複数の銘柄を、業種または企業規模毎にグループ化し、グループ毎に上記の制約条件を設定する。例えば、東京証券取引所は、約２０００銘柄を取り扱っている。また、東京証券取引所は、これら約２０００銘柄を、３３種類の業種および３種類の規模に分類している。従って、東京証券所で取り扱われている約２０００銘柄は、３３業種×３規模の合計９９個のグループに分けることができる。

【0006】

１つの銘柄の投資割合を１％から１００％までの１００段階で離散化して、東京証券取引所により取り扱われる２０００銘柄に対してポートフォリオ最適化問題を解く場合、目的関数に含まれる決定変数の個数は、２０００００個となる。

【0007】

これに対して、９９個のグループの全てに対して一次不等式による制約条件を設定するとする。上限値（Ｗ_ｍ）を１０００段階で離散化して、スラッグ変数を利用した定式化をした場合、スラッグ変数の個数は、９９０００個となる。このため、９９個のグループの全てに対して一次不等式による制約条件を設定したポートフォリオ最適化問題を解く場合、イジングスピンの個数が２００００＋９９０００＝２９９０００個のイジング問題を解かなければならない。

【0008】

このように、スラック変数を利用した定式化をしてポートフォリオ最適化問題を解く場合、制約条件が無い場合と比較して、例えば東京証券取引所の例では、イジングスピンの数が５０％程度増加してしまう。ソルバー等によりイジング問題を解く場合、イジングスピンの数が増加してしまうと、演算量および演算時間が多くなり、求解のためのコストが大きくなってしまう。

【先行技術文献】

【特許文献】

【0009】

【文献】特開２０１７－７３１０６号公報

【文献】特開２０１９－１４５０１０号公報

【文献】特開２０２１－０４３５８９号公報

【文献】国際公開第２０２０／１９６９１５号

【非特許文献】

【0010】

【文献】H. Goto,“Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network”, Sci. Rep. 6, 21686 (2016).

【文献】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

【文献】土屋和雄，西山岳宏，辻田勝吉， “分岐特性を用いた組合せ最適化問題の近似解法”， ［online］，［２０２１年５月１７日検索］，インターネット， <http://www.ynl.t.u-tokyo.ac.jp/project/RobotBrainCREST/publications/pdf/tsuchiya/4_01.pdf>

【文献】土屋和雄，西山岳宏，辻田勝吉， “分岐特性を用いた組合せ最適化問題の近似解法 第２報：決定論的アニーリングアルゴリズムの解析”， ［online］，［２０２１年５月１７日検索］，インターネット， <http://www.ynl.t.u-tokyo.ac.jp/project/RobotBrainCREST/publications/pdf/tsuchiya/4_02.pdf>

【文献】A. Lucas, “Ising formulations of many NP problems”, Frontiers in Physics 2, 5 (2014), DOI: 10.3389/fphy.2014.00005

【文献】Andrew Lucas, “Ising formulations of many NP problems”,Frontiers in Physics, Interdisciplinary Physics, Volume 2, Articles 5, 12 February 2014

【発明の概要】

【発明が解決しようとする課題】

【0011】

本発明が解決しようとする課題は、不等式により表される１または複数の制約条件の下で０－１組合せ最適化問題を少ない演算コストで解くことができる求解装置、求解方法およびプログラムを提供する。

【課題を解決するための手段】

【0012】

実施形態に係る求解装置は、０－１組合せ最適化問題を、前記０－１組合せ最適化問題に含まれる複数の離散変数を用いた不等式により表される１または複数の制約条件の下で解く。前記求解装置は、更新部と、出力部とを備える。前記更新部は、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎に順次に、前記第１変数および前記第２変数を交互に更新する。前記出力部は、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記０－１組合せ最適化問題の解を出力する。前記複数の要素は、前記複数の離散変数に対応する。前記第１変数および前記第２変数のそれぞれは、実数により表される。前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新する。前記単位時間毎の更新処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第２変数を前記第１変数に基づき更新する。前記単位時間毎の更新処理において、前記更新部は、前記１または複数の制約条件のそれぞれについて、前記複数の離散変数のそれぞれに対応する前記第１変数を代入した前記不等式を満たさない場合、前記複数の要素のそれぞれの前記第２変数から、前記不等式の境界から前記複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。

【図面の簡単な説明】

【0013】

【図1】図１は、シミュレーテッド分岐アルゴリズムの分岐現象を示す図である。

【図2】図２は、本実施形態に係る求解装置の機能構成を示す図である。

【図3】図３は、更新部の処理の流れの第１例を示すフローチャートである。

【図4】図４は、更新部の処理の流れの第２例を示すフローチャートである。

【図5】図５は、制約処理の第１例の処理の流れを示すフローチャートである。

【図6】図６は、制約処理の第２例の処理の流れを示すフローチャートである。

【図7】図７は、制約処理の第３例の処理の流れを示すフローチャートである。

【図8】図８は、情報処理システムの構成図である。

【図9】図９は、管理サーバの構成図である。

【図10】図１０は、記憶部に記憶されるデータを示す図である。

【図11】図１１は、計算サーバの構成図である。

【図12】図１２は、ストレージに記憶されるデータを示す図である。

【発明を実施するための形態】

【0014】

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

【数1】

【0015】

ｓ_ｉ、ｓ_ｊはスピンを表す。スピンは、＋１または－１のいずれかの値をとる２値変数である。ｓ_ｉは、ｉ番目のスピンを表す。ｓ_ｊは、ｊ番目のスピンを表す。ｉおよびｊは、１以上、Ｎ以下の整数である。Ｎは、スピンの数を表し、２以上の整数である。ｈ_ｉは、ｉ番目のスピンに作用する局所磁場を表す。Ｊは、２つのスピン間に作用する力を表す結合係数の行列である。Ｊは、対角成分が０である実対称行列である。Ｊ_ｉｊは、Ｊのｉ行ｊ列の要素を表す。つまり、Ｊ_ｉｊは、ｉ番目のスピンと、ｊ番目のスピンとの間に作用する力を表す結合係数である。

【0016】

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

【0017】

また、０または１のいずれかの値をとる離散変数（ビット）の２次関数を目的関数とする０－１組合せ最適化問題は、０－１二次計画問題と呼ばれる。離散変数（ビット）は、（１＋ｓ_ｉ）／２の演算を用いることにより、ｓ_ｉに変換される。つまり、０－１二次計画問題は、式（１）で表されるイジング問題と等価であるといえる。従って、０－１二次計画問題は、イジング問題に変換し、イジングマシンにより解を算出することが可能である。

【0018】

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

【0019】

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

【0020】

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

【0021】

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

【0022】

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

【数2】

【0023】

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

【数3】

【0024】

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

【0025】

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

【数4】

【0026】

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

【0027】

ｃは、係数である。ｃは、例えば、計算を実行する前に予め定められる定数であってもよい。また、α（ｔ）は、ｐ（ｔ）とともに増加する係数である。

【0028】

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

【0029】

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

【0030】

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

【数5】

【数6】

【0031】

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

【0032】

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

【0033】

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

【0034】

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

【0035】

（改良アルゴリズム）
発明者は、上述のシミュレーテッド分岐アルゴリズムを改良して、不等式により表される１または複数の制約条件の下で、０－１組合せ最適化問題を解く改良アルゴリズムを発明した。

【0036】

改良アルゴリズムは、式（３）のハミルトニアンに代えて、式（７）に示すハミルトニアンを用いて、式（２）の連立常微分方程式を数値的に解く。

【数7】

【0037】

Ｇ_ｍは、１または複数の制約条件のうちのｍ番目の制約条件に対応するペナルティ項を表す。１または複数の制約条件のそれぞれは、目的関数に含まれる複数の離散変数を含む不等式により表される。例えば、不等式は、複数の離散変数の一次項を含む一次不等式、または、複数の離散変数の二次項を含む二次不等式である。

【0038】

ペナルティ項は、対応する不等式を満たす場合に、ハミルトニアンにエネルギーとして０を加え、対応する不等式を満たさない場合に、ハミルトニアンに正のエネルギーを加える。より具体的には、Ｇ_ｍは、制約条件を表す不等式の境界から、Ｎ個の要素により特定される位置までの距離が大きい程、大きなエネルギーをハミルトニアンに加える関数である。なお、不等式の境界は、不等式を満たすか満たさないかの境界である。すなわち、不等式の境界は、不等式における不等号を、等号により代えた方程式により表される。

【0039】

そして、改良アルゴリズムは、シンプレクティック・オイラー法により、連立常微分方程式を式（８）または式（９）に示すような離散的な漸化式を用いて解く。

【数8】

【数9】

【0040】

式（８）および式（９）に示す漸化式は、式（５）および式（６）に示す漸化式と比較すると、第２変数（ｙ_ｉ（ｔ＋Δｔ））の算出式が異なる。具体的には、式（８）および式（９）に示す漸化式は、第２変数（ｙ_ｉ（ｔ＋Δｔ））の算出式に、１または複数の不等式のそれぞれのペナルティ項（Ｇ_ｍ）のｘ_ｉについての偏微分値の合計に、Δｔを乗算した値が減算されている点において異なる。

【0041】

１または複数の不等式のそれぞれの偏微分値は、対応する不等式の境界からＮ個の要素により特定される位置までの距離における対応する要素の成分に相当する。この偏微分値は、対応する不等式の境界からＮ個の要素により特定される位置までの距離における対応する要素の成分が大きい程、大きい値となる。

【0042】

ただし、ペナルティ項は、対応する不等式を満たす場合に、０となる。従って、１または複数の不等式のそれぞれの偏微分値は、Ｎ個の要素が不等式を満たしている場合には、０となる。

【0043】

ここで、シミュレーテッド分岐アルゴリズムは、時刻を表す変数であるｔを、単位時間（Δｔ）毎に増加させながら、複数の第１変数（ｘ_ｉ）および複数の第２変数（ｙ_ｉ）を交互に更新する。従って、シミュレーテッド分岐アルゴリズムは、単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれについて、不等式を満たすか否かを判断することができる。

【0044】

そこで、改良アルゴリズムは、単位時間（Δｔ）毎に、複数の第２変数（ｙ_ｉ）を更新した後に、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数（ｘ_ｉ）を代入した不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数を代入した不等式を満たさない場合、複数の第２変数（ｙ_ｉ）のそれぞれから、対応する補正値を減算する。つまり、改良アルゴリズムは、単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数を代入した不等式を満たさない場合、複数の第２変数（ｙ_ｉ）のそれぞれから、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値を減算する。

【0045】

これにより、改良アルゴリズムは、複数の第１変数（ｘ_ｉ）および複数の第２変数（ｙ_ｉ）を単位時間（Δｔ）毎に、１または複数の制約条件のそれぞれを満たすように、複数の第１変数（ｘ_ｉ）および複数の第２変数（ｙ_ｉ）を補正することができる。この結果、改良アルゴリズムは、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘ_ｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。

【0046】

このような改良アルゴリズムは、不等式により表される１または複数の制約条件の下で、０－１組合せ最適化問題を解く場合であっても、制約条件を設定しない場合と同一の変数の数により、解を得ることができる。従って、このような改良アルゴリズムは、スラッグ変数を設けずに、簡易に、不等式により表される１または複数の制約条件の下で０－１組合せ最適化問題を解くことができる。

【0047】

（制約条件が一次不等式であるのペナルティ項および偏微分値の第１例）
１または複数の制約条件の一つは、一次不等式により表される。例えば、改良アルゴリズムは、Ｍ個（Ｍは１以上の整数）の一次不等式の制約の下で、０－１組合せ最適化問題を解く。この場合において、１または複数の制約条件のうちのｍ番目の制約条件は、式（１０）により表される。

【数10】

【0048】

式（１０）において、ｍは、１以上、Ｍ以下の整数である。式（１０）において、ｘ_ｉは、ｉ番目の離散変数である。式（１０）において、ｘ_ｊは、ｊ番目の離散変数である。ω_ｍ，ｉは、ｉ番目の離散変数に乗算される係数である。Ｗ_ｍは、予め定められた定数である。

【0049】

制約条件が一次不等式により表される場合、ペナルティ項であるＧ_ｍは、式（１１）に示すように表すことができる。

【数11】

【0050】

Ａ_ｍは、ｍ番目の制約条件に対して予め定められた係数である。ｋは、１以上の整数である。例えば、ｋは、２であってもよい。

【0051】

Ｓ_ｍは、一次不等式における定数項（Ｗ_ｍ）を除く関数における複数の離散変数のそれぞれに、対応する第１変数（ｘ_ｉ）を代入して得られる値である。具体的には、Ｓ_ｍは、式（８）の漸化式を適用する場合には、式（１２）のように表される。また、Ｓ_ｍは、式（９）の漸化式を適用する場合には、式（１３）のように表される。

【数12】

【数13】

【0052】

この場合、偏微分値は、式（１４）に示すように表される。

【数14】

【0053】

そこで、制約条件が一次不等式で表される場合、改良アルゴリズムは、単位時間（Δｔ）毎に、複数の第２変数（ｙ_ｉ）を更新した後に、複数の離散変数のそれぞれに対応する第１変数（ｘ_ｉ）を代入した一次不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間（Δｔ）毎に、対応する一次不等式を満たさない場合、複数の第２変数（ｙ_ｉ）のそれぞれから、式（１４）に表される偏微分値に単位時間を表すΔｔを乗算した補正値を、減算する。

【0054】

このような処理を実行する改良アルゴリズムは、制約条件を設定しない場合と同一の変数の数により、一次不等式により表される制約条件の下で０－１組合せ最適化問題を簡易に解くことができる。

【0055】

（制約条件が一次不等式であるのペナルティ項および偏微分値の第２例）
式（１４）に示す偏微分値は、Ａ_ｍが大きい場合、計算誤差が蓄積して不安定となる可能性がある。そこで、改良アルゴリズムは、制約条件が一次不等式である場合、式（１５）により表される偏微分値を用いてもよい。すなわち、改良アルゴリズムは、単位時間（Δｔ）毎に、対応する一次不等式を満たさない場合、複数の第２変数（ｙ_ｉ）のそれぞれから、式（１５）に表される偏微分値に単位時間を表すΔｔを乗算した補正値を、減算してもよい。

【数15】

【0056】

Ｂ_ｍは、０以上１以下の値であり、ｍ番目の制約条件に対して予め定められた係数である。Ｅ_ｍは、一次不等式における定数項（Ｗ_ｍ）を除く関数における複数の離散変数のそれぞれに、対応する第２変数（ｙ_ｉ）を代入して得られる値である。具体的には、Ｅ_ｍは、式（８）または式（９）の漸化式を適用する場合には、式（１６）のように表される。

【数16】

【0057】

Ｅ_ｍは、複数の要素により特定される方向の、一次不等式の境界に対する傾きを表す。従って、式（１５）の偏微分値を用いた場合、補正値は、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分と、複数の要素により特定される方向の、一次不等式の境界に対する傾きの対応する要素の成分とを加算した値に応じた値となる。これにより、改良アルゴリズムは、Ａ_ｍが大きいために距離に応じた成分が不安定となる場合であっても、補正値を安定化させることができる。

【0058】

（制約条件が二次不等式であるのペナルティ項および偏微分値）
１または複数の制約条件の一つは、二次不等式により表されてもよい。この場合において、制約条件は、式（１７）により表される。

【数17】

【0059】

式（１７）において、ｘ_ｉは、ｉ番目の離散変数である。式（１７）において、ｘ_ｊは、ｊ番目の離散変数である。Ｑ_ｉ，ｊは、Ｎ行、Ｎ列の半正定値行列に含まれるｉ行、ｊ列の値である。ｑは、予め定められた定数である。

【0060】

制約条件が二次不等式により表される場合、ペナルティ項であるＧは、式（１８）に示すように表すことができる。

【数18】

【0061】

Ａは、二次不等式による制約条件に対して予め定められた係数である。ｋは、１以上の整数である。例えば、ｋは、２であってもよい。

【0062】

Ｓは、二次不等式における定数項（ｑ）を除く関数における複数の離散変数のそれぞれに、対応する第１変数（ｘ_ｉ，ｘ_ｊ）を代入して得られる値である。具体的には、Ｓは、式（８）の漸化式を適用する場合には、式（１９）のように表される。また、Ｓは、式（９）の漸化式を適用する場合には、式（２０）のように表される。

【数19】

【数20】

【0063】

この場合、偏微分値は、式（８）の漸化式を適用する場合には、式（２１）に示すように表される。また、偏微分値は、式（９）の漸化式を適用する場合には、式（２２）に示すように表される。

【数21】

【数22】

【0064】

そこで、制約条件が二次不等式で表される場合、改良アルゴリズムは、単位時間（Δｔ）毎に、複数の第２変数（ｙ_ｉ）を更新した後に、複数の離散変数のそれぞれに対応する第１変数（ｘ_ｉ，ｘ_ｊ）を代入した二次不等式を満たすか否かを判断する。そして、改良アルゴリズムは、単位時間（Δｔ）毎に、対応する二次不等式を満たさない場合、複数の第２変数（ｙ_ｉ）のそれぞれから、式（２１）または式（２２）に表される偏微分値に単位時間を表すΔｔを乗算した補正値を、減算する。

【0065】

イジングモデルは、二次不等式による制約条件の下で０－１組合せ最適化問題を解くことは困難である。これに対して、改良アルゴリズムは、二次不等式により表される制約条件の下で、０－１組合せ最適化問題を、簡易に解を得ることができる。

【0066】

なお、改良アルゴリズムは、１または複数の制約条件の中に、一次不等式による制約条件と、二次不等式による制約条件との両者が含まれていても、０－１組合せ最適化問題を解くことができる。また、改良アルゴリズムは、１または複数の制約条件に一次不等式による制約条件のみが含まれていてもよい。また、改良アルゴリズムは、１つの二次不等式による制約条件のみの下で、０－１組合せ最適化問題を解いてもよい。

【0067】

（機能ブロック構成）
図２は、本実施形態に係る求解装置１０の機能構成を示す図である。

【0068】

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

【0069】

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

【0070】

入力部１２は、０－１組合せ最適化問題の目的関数を定義するための情報（例えば、Ｎ、Ｊ、ｈ）、および、改良アルゴリズムを実行するために必要な係数を表す情報（例えば、Ｄ、ｃ、Δｔ、Ｔ，ｐ（ｔ）、α（ｔ））を外部装置から受け取る。さらに、入力部１２は、１または複数の制約条件を定義するための情報（例えば、ω_ｍ，ｉ，Ｗ_ｍ，Ｑ_ｉ，ｊ，ｑ）、および、制約処理を実行するために必要な係数（例えば、Ａ_ｍ，ｋ，Ｂ_ｍ，Ａ）を表す情報を外部装置から受け取る。そして、入力部１２は、受け取ったこれらの情報を更新部１４に与える。

【0071】

更新部１４は、改良アルゴリズムを用いて、第１変数（ｘ_ｉ）および第２変数（ｙ_ｉ）が対応付けられた複数の要素のそれぞれについて、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎に順次に、第１変数（ｘ_ｉ）および第２変数（ｙ_ｉ）を交互に更新する。

【0072】

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

【0073】

ここで、複数の要素は、０－１組合せ最適化問題の複数の離散変数に対応する。また、第１変数（ｘ_ｉ）および第２変数（ｙ_ｉ）のそれぞれは、実数により表される。

【0074】

そして、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第１変数（ｘ_ｉ）を第２変数（ｙ_ｉ）に基づき更新する。また、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第２変数（ｙ_ｉ）を第１変数（ｘ_ｉ）に基づき更新する。

【0075】

例えば、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第１変数（ｘ_ｉ）を更新した後に第２変数（ｙ_ｉ）を更新する。これに代えて、単位時間毎の更新処理において、更新部１４は、複数の要素のそれぞれについて、第２変数（ｙ_ｉ）を更新した後に第１変数（ｘ_ｉ）を更新してもよい。

【0076】

さらに、単位時間毎の更新処理において、更新部１４は、第２変数（ｙ_ｉ）を更新した後に、０－１組合せ最適化問題の解が１または複数の制約条件を満たすように、制約処理を実行する。更新部１４は、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数（ｘ_ｉ）を代入した不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙ_ｉ）から、対応する不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値（ｒ_ｍ，ｉ）を減算する。

【0077】

単位時間毎の更新処理において、１または複数の制約条件のそれぞれについて制約処理を実行することにより、更新部１４は、複数の第１変数（ｘ_ｉ）が制約条件を満たしていない場合、複数の第１変数（ｘ_ｉ）が制約条件を満たす方向に変化するように、複数の要素に対応する複数の第２変数（ｙ_ｉ）を補正することができる。これにより、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘ_ｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。

【0078】

なお、単位時間毎の更新処理において、更新部１４は、第１変数（ｘ_ｉ）を更新した後に第２変数（ｙ_ｉ）を更新する場合には、複数の要素のそれぞれについて、第２変数（ｙ_ｉ）を更新した後に第２変数（ｙ_ｉ）から補正値（ｒ_ｍ，ｉ）を減算する。また、単位時間毎の更新処理において、更新部１４は、第２変数（ｙ_ｉ）を更新した後に第１変数（ｘ_ｉ）を更新する場合には、複数の要素のそれぞれについて、第２変数（ｙ_ｉ）を更新した後、且つ、第１変数（ｘ_ｉ）を更新する前に、第２変数（ｙ_ｉ）から補正値（ｒ_ｍ，ｉ）を減算する。

【0079】

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

【0080】

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

【0081】

続いて、Ｓ１０２において、更新部１４は、制約条件に関する情報を設定する。具体的には、更新部１４は、１または複数の制約条件を定義するための情報、および、制約処理を実行するために必要な係数を表す情報を設定する。例えば、更新部１４は、制約条件が一次不等式により表される場合、変数に乗算される係数であるω_ｍ，ｉ、定数項であるＷ_ｍ、制約処理を実行するために必要な係数であるＡ_ｍ（またはＡ_ｍおよびＢ_ｍ）、並びに、ｋを設定する。また、例えば、制約条件が二次不等式により表される場合、更新部１４は、変数に乗算される係数であるＱ_ｉ，ｊおよび定数項であるｑ、制約処理を実行するために必要な係数であるＡおよびｋを設定する。なお、更新部１４は、制約処理を実行するために必要な係数であるＡ_ｍ、Ｂ_ｍ、ｋおよびＡを、入力部１２から受け取ったパラメータに応じて設定してもよいし、予め決定されており変更できない値を設定してもよい。

【0082】

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

【0083】

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

【0084】

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

【0085】

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

【0086】

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

【数23】

【0087】

すなわち、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））と、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））とに基づき更新する。

【0088】

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

【0089】

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

【0090】

Ｓ１０９において、更新部１４は、Ｎ個の要素のそれぞれの対象時刻（ｔ＋Δｔ）における第１変数（ｘ_１（ｔ＋Δｔ）～ｘ_Ｎ（ｔ＋Δｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に０－１組合せ最適化問題により予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ＋Δｔ））を算出する。作用係数は、Ｊに含まれる結合係数およびｈに含まれる局所磁場係数である。具体的には、更新部１４は、式（２４）を算出する。

【数24】

【0091】

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

【数25】

【0092】

続いて、Ｓ１１１において、更新部１４は、時間経過に従って増加する関数であるｐ（ｔ＋Δｔ）に基づき定まる値に、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_１（ｔ＋Δｔ））を乗算した時間発展値（ｇ_ｉ（ｔ＋Δｔ））を算出する。具体的には、更新部１４は、式（２６）を算出する。

【数26】

【0093】

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

【数27】

【0094】

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

【0095】

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

【0096】

Ｓ１１４において、更新部１４は、１または複数の制約条件に基づく制約処理を実行する。なお、制約処理については、図５～図７のフローチャートを用いて詳細を後述する。Ｓ１１４の処理を終えると、更新部１４は、処理をＳ１１５に進める。

【0097】

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

【0098】

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

【0099】

以上のＳ１０１～Ｓ１１６の処理を実行することにより、更新部１４は、改良アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

【0100】

図４は、更新部１４の処理の流れの第２例を示すフローチャートである。更新部１４は、改良アルゴリズムを用いて０－１組合せ最適化問題を解く場合、図３に示す流れに代えて、図４に示す流れで処理を実行してもよい。

【0101】

まず、Ｓ２０１、Ｓ２０２およびＳ２０３において、更新部１４は、図３に示す第１例のＳ１０１、Ｓ１０２およびＳ１０３と同一の処理を実行する。

【0102】

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

【0103】

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

【0104】

Ｓ２０６において、更新部１４は、Ｎ個の要素のそれぞれの直前時刻（ｔ）における第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に０－１組合せ最適化問題により予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ））を算出する。具体的には、更新部１４は、式（２８）を算出する。

【数28】

【0105】

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

【数29】

【0106】

続いて、Ｓ２０８において、更新部１４は、時間経過に従って増加する関数であるｐ（ｔ）に基づき定まる値に、対象要素の直前時刻（ｔ）における第１変数（ｘ_１（ｔ））を乗算した時間発展値（ｇ_ｉ（ｔ））を算出する。具体的には、更新部１４は、式（３０）を算出する。

【数30】

【0107】

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

【数31】

【0108】

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

【0109】

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

【0110】

Ｓ２１１において、更新部１４は、１または複数の制約条件に基づく制約処理を実行する。なお、制約処理については、図５～図７のフローチャートを用いて詳細を後述する。Ｓ２１１の処理を終えると、更新部１４は、処理をＳ２１２に進める。

【0111】

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

【0112】

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

【数32】

【0113】

すなわち、更新部１４は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））と、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））とに基づき更新する。

【0114】

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

【0115】

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

【0116】

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

【0117】

以上のＳ２０１～Ｓ２１６の処理を実行することにより、更新部１４は、改良アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

【0118】

（制約処理の第１例）
図５は、制約処理の第１例の処理の流れを示すフローチャートである。制約条件として、１または複数の制約条件としてＭ個の一次不等式が設定されている場合、更新部１４は、図３のＳ１１４および図４のＳ２１１において、図５に示す流れで処理を実行する。

【0119】

まず、更新部１４は、Ｓ３０１とＳ３０８との間のループ処理を、ｍ＝１からｉ＝Ｍまでｍを１ずつインクリメントしながら繰り返す。なお、ｍは、Ｍ個の一次不等式のうちの処理対象となる一次不等式を表すインデックスである。

【0120】

Ｓ３０１とＳ３０８とのループ内において、まず、Ｓ３０２において、更新部１４は、制約値（Ｓ_ｍ）を算出する。具体的には、図３のＳ１１４の制約処理において、更新部１４は、式（３３）の演算を実行する。また、図４のＳ２１１の制約処理において、更新部１４は、式（３４）の演算を実行する。

【数33】

【数34】

【0121】

なお、制御値であるＳ_ｍは、対応する一次不等式における定数項（Ｗ_ｍ）を除く関数に含まれる複数の離散変数のそれぞれに、対応する第１変数（ｘ_ｉ）を代入し得られる値である。

【0122】

続いて、Ｓ３０３において、更新部１４は、制約値（Ｓ_ｍ）が対応する一次不等式における定数項（Ｗ_ｍ）以下であるか否かを判断する。すなわち、更新部１４は、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たすか否かを判断する。

【0123】

更新部１４は、制約値（Ｓ_ｍ）が対応する一次不等式における定数項（Ｗ_ｍ）以下である場合、すなわち、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たす場合（Ｓ３０３のＹｅｓ）、処理をＳ３０８に進める。更新部１４は、制約値（Ｓ_ｍ）が対応する一次不等式における定数項（Ｗ_ｍ）以下ではない場合、すなわち、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たさない場合（Ｓ３０３のＮｏ）、処理をＳ３０４に進める。

【0124】

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

【0125】

Ｓ３０４とＳ３０７とのループ内において、まず、Ｓ３０５において、更新部１４は、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値（ｒ_ｍ，ｉ）を算出する。具体的には、更新部１４は、式（３５）の演算を実行する。

【数35】

【0126】

続いて、Ｓ３０６において、更新部１４は、対象要素の第２変数（ｙ_ｉ）から、補正値（ｒ_ｍ，ｉ）を減算する。具体的には、更新部１４は、式（３６）の演算を実行する。

【数36】

【0127】

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

【0128】

そして、更新部１４は、Ｓ３０１とＳ３０８との間のループ処理をＭ回実行した場合、すなわち、Ｍ個の一次不等式の全てについて処理を実行した場合、本フローを終了する。

【0129】

以上の処理を実行することにより、更新部１４は、Ｍ個の一次不等式のそれぞれについて、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙ_ｉ）から、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値（ｒ_ｍ，ｉ）を減算することができる。これにより、更新部１４は、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たしていない場合、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たす方向に変化するように、複数の第２変数（ｙ_ｉ）を補正することができる。この結果、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘ_ｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。

【0130】

（制約処理の第２例）
図６は、制約処理の第２例の処理の流れを示すフローチャートである。制約条件として、１または複数の制約条件としてＭ個の一次不等式が設定されている場合、更新部１４は、図３のＳ１１４および図４のＳ２１１において、図６に示す流れで処理を実行してもよい。

【0131】

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

【0132】

Ｓ４０１とＳ４０９とのループ内において、まず、Ｓ４０２において、更新部１４は、制約値（Ｓ_ｍ）を算出する。Ｓ４０２の処理は、図５のＳ３０２の処理と同一である。

【0133】

続いて、Ｓ４０３において、更新部１４は、制約値（Ｓ_ｍ）が対応する一次不等式における定数項（Ｗ_ｍ）以下であるか否かを判断する。Ｓ４０３の処理は、図５のＳ３０３の処理と同一である。更新部１４は、制約値（Ｓ_ｍ）が対応する一次不等式における定数項（Ｗ_ｍ）以下である場合（Ｓ４０３のＹｅｓ）、処理をＳ４０９に進める。更新部１４は、制約値（Ｓ_ｍ）が対応する一次不等式における定数項（Ｗ_ｍ）以下ではない場合（Ｓ４０３のＮｏ）、処理をＳ４０４に進める。

【0134】

Ｓ４０４において、更新部１４は、複数の要素により特定される方向における対応する一次不等式の境界に対する傾き（Ｅ_ｍ）を算出する。具体的には、更新部１４は、式（３７）の演算を実行する。

【数37】

【0135】

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

【0136】

Ｓ４０５とＳ４０８とのループ内において、まず、Ｓ４０６において、更新部１４は、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分、および、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた補正値（ｒ_ｍ，ｉ）を算出する。具体的には、更新部１４は、式（３８）の演算を実行する。

【数38】

【0137】

続いて、Ｓ４０７において、更新部１４は、対象要素の第２変数（ｙ_ｉ）から、補正値（ｒ_ｍ，ｉ）を減算する。Ｓ４０７の処理は、図５のＳ３０６の処理と同一である。

【0138】

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

【0139】

そして、更新部１４は、Ｓ４０１とＳ４０９との間のループ処理をＭ回実行した場合、すなわち、Ｍ個の一次不等式の全てについて処理を実行した場合、本フローを終了する。

【0140】

以上の処理を実行することにより、更新部１４は、Ｍ個の一次不等式のそれぞれについて、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙ_ｉ）から、一次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分、および、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた補正値（ｒ_ｍ，ｉ）を減算することができる。特に、更新部１４は、補正値（ｒ_ｍ，ｉ）に、複数の要素により特定される方向の一次不等式の境界に対する傾きの対象要素の成分に応じた値を加えることにより、単位時間（Δｔ）毎の更新処理において、複数の第１変数（ｘ_ｉ）および複数の第２変数（ｙ_ｉ）が不安定に変動することを抑制することができる。

【0141】

これにより、更新部１４は、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たしていない場合、複数の第１変数（ｘ_ｉ）が対応する一次不等式を満たす方向に安定して変化するように、複数の要素に対応する複数の第２変数（ｙ_ｉ）を補正することができる。この結果、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘ_ｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。

【0142】

（制約処理の第３例）
図７は、制約処理の第３例の処理の流れを示すフローチャートである。１または複数の制約条件のうちの一つに二次不等式が設定されている場合、更新部１４は、二次不等式については、図５のＳ３０２～Ｓ３０７または図６のＳ４０２～Ｓ４０８の処理に代えて、図７に示すＳ５０１～Ｓ５０６の処理を実行する。

【0143】

まず、Ｓ５０１において、更新部１４は、制約値（Ｓ）を算出する。具体的には、図３のＳ１１４の制約処理において、更新部１４は、式（３９）の演算を実行する。また、図４のＳ２１１の制約処理において、更新部１４は、式（４０）の演算を実行する。

【数39】

【数40】

【0144】

なお、制御値であるＳは、二次不等式における定数項（ｑ）を除く関数に含まれる複数の離散変数のそれぞれに、対応する第１変数（ｘ_ｉ）を代入して得られる値である。

【0145】

続いて、Ｓ５０２において、更新部１４は、制約値（Ｓ）が二次不等式における定数項（ｑ）以下であるか否かを判断する。すなわち、更新部１４は、複数の第１変数（ｘ_ｉ）が二次不等式を満たすか否かを判断する。

【0146】

更新部１４は、制約値（Ｓ）が二次不等式における定数項（ｑ）以下である場合、すなわち、複数の第１変数（ｘ_ｉ）が対応する二次不等式を満たす場合（Ｓ５０２のＹｅｓ）、本フローを終了する。更新部１４は、制約値（Ｓ）が二次不等式における定数項（ｑ）以下ではない場合、すなわち、複数の第１変数（ｘ_ｉ）が対応する二次不等式を満たさない場合（Ｓ５０２のＮｏ）、処理をＳ５０３に進める。

【0147】

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

【0148】

Ｓ５０３とＳ５０６とのループ内において、まず、Ｓ５０４において、更新部１４は、二次不等式の境界から複数の要素により特定される位置までの距離における対象要素の成分に応じた補正値（ｒ_ｉ）を算出する。具体的には、図３のＳ１１４の制約処理において、更新部１４は、式（４１）の演算を実行する。また、図４のＳ２１１の制約処理において、更新部１４は、式（４２）の演算を実行する。

【数41】

【数42】

【0149】

続いて、Ｓ５０５において、更新部１４は、対象要素の第２変数（ｙ_ｉ）から、補正値（ｒ_ｉ）を減算する。具体的には、更新部１４は、式（４３）の演算を実行する。

【数43】

【0150】

更新部１４は、Ｓ５０３とＳ５０６との間のループ処理をＮ回実行した場合、本フローを終了する。

【0151】

以上の処理を実行することにより、更新部１４は、二次不等式について、複数の第１変数（ｘ_ｉ）が対応する二次不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙ_ｉ）から、二次不等式の境界から複数の要素により特定される位置までの距離の対象要素の成分に応じた補正値（ｒ_ｉ）を減算することができる。これにより、更新部１４は、複数の第１変数（ｘ_ｉ）が対応する二次不等式を満たしていない場合、複数の第１変数（ｘ_ｉ）が対応する二次不等式を満たす方向に変化するように、複数の第２変数（ｙ_ｉ）を補正することができる。この結果、更新部１４は、終了時刻（ｔ＝Ｔ）における複数の第１変数（ｘ_ｉ）が１または複数の制約条件の全てを満たす確率を高くすることができる。

【0152】

以上のように、本実施形態に係る求解装置１０は、シミュレーテッド分岐アルゴリズムにより０－１組合せ最適化問題を解く。本実施形態に係る求解装置１０は、単位時間（Δｔ）毎の更新処理において、１または複数の制約条件のそれぞれについて、複数の離散変数のそれぞれに対応する第１変数を代入した不等式を満たさない場合、複数の要素のそれぞれの第２変数（ｙ_ｉ）から、不等式の境界から複数の要素により特定される位置までの距離における対応する要素の成分に応じた補正値（ｒ_ｉ）を減算する。

【0153】

これにより、本実施形態に係る求解装置１０は、不等式により表される１または複数の制約条件の下で、０－１組合せ最適化問題を解く場合であっても、制約条件を設定しない場合と同一の変数の数により、解を得ることができる。従って、本実施形態に係る求解装置１０によれば、スラッグ変数を設けずに、簡易に、不等式により表される１または複数の制約条件の下で０－１組合せ最適化問題を解くことができる。この結果、本実施形態に係る求解装置１０は、変数の途中経過等を記憶するメモリ量を少なくし、計算量および計算時間を少なくすることができる。

【0154】

なお、求解装置１０は、一次不等式の制約処理を実行する場合、Ｍ個の一次不等式を生成し、ｉ番目の離散変数に乗算される係数であるω_ｍ，ｉ、および、一次不等式の定数であるＷ_ｍを呼び出してもよい。

【数44】

【0155】

求解装置１０は、一回の列番号および行番号の指定により、ω_ｍ，ｉおよびＷ_ｍを同時に読み出すことができ、容易に処理を実行することができる。

【0156】

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

【0157】

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

【0158】

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

【0159】

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

【0160】

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

【0161】

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

【0162】

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

【0163】

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

【0164】

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

【0165】

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

【0166】

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

【0167】

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

【0168】

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

【0169】

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

【0170】

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

【0171】

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

【0172】

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

【0173】

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

【0174】

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

【0175】

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

【0176】

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

【0177】

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

【0178】

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

【符号の説明】

【0179】

１０ 求解装置
１２ 入力部
１４ 更新部
１６ 出力部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】