JP2022052222A - 最適化装置、最適化方法、および最適化プログラム - Google Patents

Abstract

【課題】複数のレプリカを用いた場合の解探索能力を向上させる。【解決手段】最適化装置は、複数の状態変数を有する複数のレプリカそれぞれについて、レプリカが有する複数の状態変数のうちの第１状態変数の値を更新した場合における、そのレプリカと、複数のレプリカからそのレプリカを除いたレプリカ群のうちの一部である他のレプリカとの、複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での距離の変化に応じた相互作用の強さの変化量を特定する。そして最適化装置は、第１状態変数の値を更新した場合における相互作用の強さの変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、第１状態変数の値を更新するか否かを決定する。【選択図】図１９

Description

本発明は、最適化装置、最適化方法、および最適化プログラムに関する。

ノイマン型コンピュータが不得意とする問題として、大規模な離散最適化問題がある。離散最適化問題を計算する装置としては、例えば、イジング型の評価関数（エネルギー関数などとも呼ばれる）を用いたイジングマシン（ボルツマンマシンとも呼ばれる）がある（例えば、特許文献１参照）。

イジングマシンによる計算では、計算対象の問題は磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えられる。そして、マルコフ連鎖モンテカルロ法により、イジング型の評価関数の値（イジングモデルのエネルギーに相当する）が最小となる状態の探索が行われる。以下、マルコフ連鎖モンテカルロ法を、ＭＣＭＣ（Markov-Chain Monte Carlo）法と略す。ＭＣＭＣ法では、例えばメトロポリス法またはギブス法で規定される状態遷移の受け入れ確率で、その状態遷移が受け入れられる。

ＭＣＭＣ法の一種として、レプリカ交換法（交換モンテカルロ法またはパラレルテンパリング（parallel tempering）法とも呼ばれる）がある。レプリカ交換法は複数の温度を用いたＭＣＭＣ処理を互いに独立に行い、ある試行回数ごとに、各ＭＣＭＣ処理で得られるエネルギーを比較し、適切な確率で２つの温度に対する状態を交換するという操作を行う方法である。レプリカ交換によれば、温度を徐々に下げていく疑似焼き鈍し法と比べて、局所解に拘束される可能性が抑えられ、全探索空間を効率よく探索できる。

なお、従来、回路の物量を削減しつつ、メトロポリス法に基づく確率的な処理を可能とする情報処理装置が提案されている（例えば、特許文献１参照）。また、分子動力学シミュレーションの分野において、２つの分子間の位相距離を計算して分子間の相互作用を除外するか否かを抑制させるか否かを決定する技術が提案されている（例えば、特許文献２参照）。複数のレプリカを用いて解探索を行う技術としては、Collective Monte Carlo（ＣＭＣ）と呼ばれる方法や、Robust Ensemble（ＲＥ）と呼ばれる方法も提案されている（例えば、非特許文献１，２参照）。

特開２０１９－０８２７９３号公報 米国特許出願公開第２０１９／００８７５４６号明細書

Gregoire Clarte and Antoine Diez,"Collective sampling through a Metropolis-Hastings like method: kinetic theory and numerical experiments", arXiv:1909.08988v1 [math.ST], 18 Sep. 2019 Baldassi, Carlo. et. al., "Unreasonable Effectiveness of Learning Neural Networks: From Accessible States and Robust Ensembles to Basic Algorithmic Schemes", PNAS E7655-E7662, Published online 15 Nov. 2016

ＭＣＭＣ法を高速化するために、多数のレプリカを用いて集団で探索をする様々な方法が提案されているが、いずれの方法においても、集団探索の効果が十分には発揮できない場合がある。例えば、遷移先の候補の選び方が１－ｂｉｔフリップ（複数のビットのうちの１つの値を反転させる）の場合、各ビットが反転対象として等確率で選択され、選択されたビットを反転させた状態への遷移確率は遷移前後のエネルギー差に基づいて決定される。そのため、どのレプリカもエネルギー勾配に従って状態変化し、状態遷移の過程が同じ道筋を辿ってしまう可能性がある。その結果、複数のレプリカが同じ局所解に留まってしまい、十分に広く状態空間を探索できない場合がある。

なお、このような問題は、状態変数の値が離散的な場合に限らず、状態変数が連続的な値を取り得る最適化問題においても同様に生じる。
１つの側面では、本発明は、複数のレプリカを用いた場合の解探索能力を向上させることを目的とする。

１つの案では、以下に示す記憶部と処理部とを有する最適化装置が提供される。
記憶部は、複数のレプリカそれぞれの複数の状態変数の値を記憶する。処理部は、複数のレプリカそれぞれについて、レプリカが有する複数の状態変数のうちの第１状態変数の値を更新した場合における、そのレプリカと、複数のレプリカからそのレプリカを除いたレプリカ群のうちの一部である他のレプリカとの、複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での距離の変化に応じた相互作用の強さの変化量を特定する。そして処理部は、第１状態変数の値を更新した場合における相互作用の強さの変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、第１状態変数の値を更新するか否かを決定する。

１態様によれば、複数のレプリカを用いた場合の解探索能力を向上させることができる。

本実施の形態に係る最適化方法の比較例を示す図である。 第１の実施の形態に係る最適化方法の一例を示す図である。 第２の実施の形態のシステム構成の一例を示す図である。 サーバのハードウェアの一例を示す図である。 イジングマシンの一例を示す図である。 イジングモデルの模式図である。 レプリカ交換の一例を示す図である。 Ｏｎｅ－Ｈｏｔ制約下での１－ｂｉｔフリップの例を示す図である。 ２Ｗ１Ｈの制約を説明する図である。 イジングマシンの解探索機能の一例を示す図である。 解探索エンジンにおける処理の一例を示す図である。 相互作用を与えるレプリカの１つ目の選択方法の一例を示す図である。 相互作用を与えるレプリカの１つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。 レプリカごとの解探索処理の手順の一例を示すフローチャートである。 レプリカ間の相互作用のエネルギーの差分の計算手順の一例を示すフローチャートである。 相互作用を与えるレプリカの２つ目の選択方法の一例を示す図である。 相互作用を与えるレプリカの２つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。 相互作用を与えるレプリカの３つ目の選択方法の一例を示す図である。 相互作用を与えるレプリカの３つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。 相互作用を与えるレプリカの３つ目の選択方法におけるレプリカごとの解探索処理の手順の一例を示すフローチャートである。 レプリカ間の相互作用のエネルギーの差分の計算手順の一例を示すフローチャートである。 相互作用を与えるレプリカの４つ目の選択方法の一例を示す図である。 相互作用を与えるレプリカの４つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。 相互作用を与えるレプリカの４つ目の選択方法におけるレプリカごとの解探索処理の手順の一例を示すフローチャートである。 第１の更新ビット選択方法による更新ビット選択処理の手順の一例を示すフローチャートである。 第２の更新ビット選択方法の処理手順の一例を示すフローチャートである。 更新ビット選択のためのツリー状に接続されたセレクタの一例を示す図である。 第３の更新ビット選択方法の処理手順の一例を示すフローチャートである。 レプリカ間に斥力の相互作用を設定した場合のエネルギーランドスケープを示す図である。 レプリカ間に引力の相互作用を設定した場合のエネルギーランドスケープを示す図である。 第１の検証例を示す図である。 第２の検証例を示す図である（その１）。 第２の検証例を示す図である（その２）。

以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔比較例〕
まず、本実施の形態に係る最適化方法の比較例を説明する。

図１は、本実施の形態に係る最適化方法の比較例を示す図である。
図１には、解探索方法を実現する最適化装置１０ａが示されている。最適化装置１０ａは、ノイマン型コンピュータであってもよく、非ノイマン型コンピュータであってもよい。例えば最適化装置１０ａは、最適化用のプログラムを実行することにより、最適化方法を実施することができる。また最適化装置１０ａは、イジングモデルを用いた最適化問題の求解を行うイジングマシンであってもよい。イジングマシンには、量子ビットを用いた量子コンピュータ、量子ビットの量子現象をデジタル回路上で再現した装置などが含まれる。

最適化装置１０ａは、記憶部１１ａと処理部１２ａとを有する。記憶部１１ａは、例えば最適化装置１０ａが有するメモリ、またはストレージ装置である。処理部１２ａは、例えば最適化装置１０ａが有するプロセッサ、または演算回路である。演算回路には、量子ビットまたは量子ビットの仕組みを再現するニューロン回路が含まれる。

記憶部１１ａは、複数のレプリカ２～４それぞれの複数の状態変数の値を記憶する。
処理部１２ａは、複数のレプリカ２～４を用いて、最適化問題を求解する。例えば処理部１２ａは、最適化問題に応じて定義された目的関数の値が最小となる状態変数の値を求める。目的関数は、最適化問題を表すモデルのエネルギーと呼ばれることもある。最適化問題がイジングモデルで表される場合、イジングモデルのハミルトニアンが、エネルギーを示す目的関数に相当する。

解探索のために処理部１２ａは、複数のレプリカ２～４それぞれについて状態遷移（状態変数の値を更新すること）を繰り返し、生成された状態における複数の状態変数の値に基づいて目的関数の値を計算する。その際、処理部１２ａは、レプリカ間の相互作用を考慮して、レプリカの状態遷移を行う。レプリカ間の相互作用としては、例えばレプリカ間の距離に応じた引力または斥力が考えられる。ｋ番目のレプリカｘ^ｋとｌ番目のレプリカｘ^ｌとの距離は、ｄ（ｘ^ｋ，ｘ^ｌ）と表記する（ｋ，ｌは、１以上の整数）。例えば処理部１２ａは、複数のレプリカ２～４それぞれについて以下のようにして状態遷移を行う。

処理部１２ａは、レプリカが有する複数の状態変数のうちの第１状態変数の値を更新した場合における、複数の状態変数の値の組み合わせが取り得る空間を示す状態空間１内でのそのレプリカと他のレプリカとの距離の変化に応じた相互作用の強さの変化量を特定する。相互作用の強さは、例えば他のレプリカそれぞれとの距離の合計に基づく値である。相互作用の強さは、相互作用のエネルギーＧ（Ｘ）と呼ぶこともできる。相互作用の強さは、例えば後述の式（１５）または式（１６）で表すことができる。ｌ番目のレプリカのｊ_０番目の状態変数を更新した場合の相互作用の強さの変化量は、ΔＧ＝Ｇ（ｘ^ｌ［ｊ_０］）－Ｇ（ｘ^ｌ）と表すことができる。

そして処理部１２ａは、第１状態変数（例えばｊ_０番目の状態変数）の値を更新した場合における、第１状態変数の値を更新するか否かを決定する。この決定は、相互作用の強さの変化量に応じた提案確率（ｇ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）と、目的の確率分布に応じた受け入れ確率（ａ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）とに基づいて、確率的に行われる。目的の確率分布は、例えばギブス分布である。提案確率と受け入れ確率とに基づいてレプリカが状態遷移する遷移確率は、例えばメトロポリスヘイスティング法に従う。

処理部１２ａは、第１状態変数の値を更新すると決定した場合、第１状態変数の値を更新した後のレプリカの複数の状態変数の値に基づいて、目的関数の値を計算する。また処理部１２ａは、記憶部１１ａ内におけるレプリカの第１状態変数の値を更新する。そして処理部１２ａは、複数のレプリカ２～４それぞれの複数の状態変数の一の状態変数の値の更新を繰り返し、目的関数の値が所定の条件を満たしたときの複数の状態変数の値を出力する。例えば処理部１２は、複数のレプリカ２～４の更新を所定回数繰り返した後、目的関数の値が最小となる複数の状態変数の値の組み合わせを出力する。

このようにして、レプリカ間の相互作用を考慮したレプリカの状態遷移により、解の探索が行われる。すなわち最適化装置１０ａは、レプリカ間の相互作用を考慮したことで、複数のレプリカ２～４による状態空間１内を網羅的に探索することができる。しかもメトロポリスヘイスティング法を用いることで、最適化装置１０ａは、レプリカ間の相互作用による影響を適切な形で計算に組み込むことができる。

なお処理部１２ａは、状態空間１に対して適切な距離を定義し、レプリカ間の距離を定める。そして処理部１２ａは、その距離を用いてレプリカ間の相互作用の強さを決定し、メトロポリスヘイスティング法における遷移候補先の分布（提案分布）を定義し、計算に組み込む。メトロポリスヘイスティング法は提案分布が非対称な場合に対応している。そのため、提案分布の決め方に自由度がある。そこで、処理部１２ａは、メトロポリスヘイスティング法における提案分布（提案確率の定義）の自由度を利用し、提案確率内にレプリカ間相互作用を導入している。

レプリカ間の相互作用として、例えば斥力の相互作用を生じさせることができる。この場合、処理部１２ａは、第１状態変数の値を更新すると状態遷移判断対象のレプリカと他のレプリカとの距離が遠ざかる場合に、相互作用の強さを増加させる。処理部１２ａは、相互作用の強さの増加量が大きいほど提案確率を大きくする。そして処理部１２ａは、提案確率が大きい状態変数ほど、値を更新する候補として選択される確率を高くする。その結果、例えば探索空間において複数のレプリカ２～４を分散させて効率的に探索したり、複数のレプリカ２～４が、同じ局所解に嵌まり出られなくなることを抑止できる。

またレプリカ間の相互作用として、引力の相互作用を生じさせることもできる。この場合、処理部１２ａは、第１状態変数の値を更新すると状態遷移判断対象のレプリカと他のレプリカとの距離が近づく場合に、相互作用の強さを増加させる。処理部１２は、相互作用の強さの増加量が大きいほど提案確率を大きくする。これにより、例えば探索空間の特定の空間を複数のレプリカ２～４を用いて集中的に探索したり、局所解に嵌まって出られなくなっているレプリカを、他のレプリカからの引力により、局所解から抜け出させることが可能となる。

状態空間１が離散的であり、状態変数の値が二値（例えば「１」または「０」）のみを取り得る場合、２つのレプリカ間の距離として、例えばハミング距離（またはその単調増加関数）を用いることができる。この場合、処理部１２ａは、すべてのレプリカ間のハミング距離を定義し、それによりレプリカ相互作用の強さを計算する。レプリカ間の距離は、後述の式（１９）のように表すことができる。

なお処理部１２ａは、第１状態変数の値を更新することについての提案確率を、例えば規格化定数によって規格化する。例えば第１状態変数の値を更新した場合における相互作用の強さの変化量をΔＧ、レプリカに設定された温度パラメータの値の逆数である逆温度をβとする。このとき処理部１２ａは、ｅｘｐ（－βΔＧ）を所定の規格化定数で除算した値を提案確率とする。この提案確率は、例えば後述の式（１７）で表すことができる。ｅｘｐ（－βΔＧ）はギブス分布を表しており、提案確率の定義にギブス分布を用いたことで、目的関数（エネルギー）においてギブス分布を保つことが容易になる。

また処理部１２ａは、１とｅｘｐ（－βΔＧ）とのうちの小さい方を、所定の規格化定数で除算した値を提案確率としてもよい。この提案確率は、例えば後述の式（１８）で表すことができる。これによりｅｘｐ（－βΔＧ）が１を超える場合は１とみなされることとなり、相互作用の強さの変化量が状態変数間で大きく異なる場合の提案確率への影響の差を弱めることができる。

ここで、規格化定数について説明する。従来の提案分布は複数の状態変数が等確率（１／Ｎ）で遷移候補として選択される（Ｎは、状態変数の数を示す１以上の整数）。この場合、規格化定数はＮ（それぞれの遷移先の重みが共通で１ということ）である。図１の最適化装置１０ａでは、遷移候補となる状態変数それぞれの遷移確率が異なり、規格化定数が遷移前の現在状態に依存する。そこで処理部１２ａにおいて規格化定数を計算することとなる。

例えば処理部１２ａは、複数の状態変数それぞれを第１状態変数とした場合における複数の状態変数それぞれについてのｅｘｐ（－βΔＧ）の値の総和を、規格化定数とする。この規格化定数は、例えば後述の式（２３）で表すことができる。なお、相互作用がハミング距離の一次関数の場合、処理部１２ａは、レプリカの状態遷移ごとに状態遷移前と状態遷移後との規格化定数の差分計算を行い、差分の累積値を計算し（累積計算）、最新の規格化定数とすることができる。ハミング距離の一次関数は、後述の式（１９）に示すような関数である。

規格化定数の累積計算をする場合、処理部１２ａは、レプリカを状態遷移させるごとに、更新させる状態変数の決定に使用した規格化定数を記憶部１１に格納しておく。そして処理部１２ａは、レプリカの状態遷移の際に使用した規格化定数の値と、前回の状態遷移の前後で生じる規格化定数の値の差分とに基づいて、今回の状態遷移で使用する規格化定数の値を算出する。前回の状態遷移の前後で生じる規格化定数の値の差分は、例えば後述の式（２４）で表される。これにより、規格化定数を効率的に算出することができる。

なお処理部１２ａは、他のレプリカそれぞれとの距離の平方根の合計に基づく値を、相互作用の強さとしてもよい。この場合の相互作用の強さは、例えば後述の式（１６）で表される。これにより、距離が遠い他のレプリカよりも距離が近い他のレプリカからの相互作用が、相対的により強く働くようにすることができる。例えば、複数のレプリカ２～４が同じ局所解に嵌まることを抑止する場合、その局所解の近辺に存在するレプリカ間に強い斥力を働かせることで、局所解からの脱出を促進することができる。この場合、その局所解とは遠く離れた位置のレプリカからの影響が少ない方が、局所解からの脱出が容易となる。

また処理部１２ａは、複数の状態変数の中から、値の更新を受け入れることができる状態変数を先に特定し、その中から今回のレプリカの状態遷移において値を更新する状態変数を決定することもできる。この場合、処理部１２ａは、複数の状態変数それぞれについて、受け入れ確率に基づいて、状態変数の更新が提案された場合に更新を受け入れるか否かを確率的に決定する。そして処理部１２ａは、更新を受け入れると判定された状態変数の中から、提案確率が高い状態変数ほど選択される可能性を高くして、少なくとも１つの状態変数を更新対象に決定する。これにより、選択した状態変数の値の更新の棄却（更新を受け入れないとの判定）が繰り返され、値を更新する状態変数の決定に時間がかかることを抑止できる。

ところで、上記の比較例では、処理部１２ａは、全てのレプリカ間の相互作用の強さを考慮して、相互作用の強さの変化量を計算する。ただ、全てのレプリカ間の相互作用を考慮しなくても、局所解からの脱出効果が得られる場合が多く、全てのレプリカ間の相互作用を考慮した場合、むしろ状態遷移を阻害する可能性もある。例えば、すでに状態が大きく異なっている（距離が遠い）レプリカ間には、斥力の相互作用を発生させる必要性が低い。状況によっては、距離が遠いレプリカの影響を大きく受けて、状態遷移が阻害されることがある。また、引力の相互作用を発生させる場合、複数のレプリカが同じ局所解に嵌まってしまうと、他のレプリカをより強くその局所解に引き付けてしまうため、状態遷移を阻害してしまうことがある。このため、相互作用を与えるレプリカの範囲を制限したほうがよい場合がある。

また、相互作用の強さを表す後述の式（１５）または式（１６）に含まれるレプリカ間の距離の計算回数は、全レプリカ数をＭ（Ｍは２以上の整数）とすると、各レプリカについてＭ回であるため、最適化装置１０ａ全体では、Ｍ^２回となる。たとえば、Ｍ＝１００の場合、上記の計算回数は１０^４回、Ｍ＝１０００の場合、計算回数は１０^６回となる。このため、レプリカ数が増えるほど計算量が大幅に増大する。

以下に示す第１の実施の形態に係る最適化方法は、上記比較例に対して、計算量を抑えることを可能とするものである。
〔第１の実施の形態〕
図２は、第１の実施の形態に係る最適化方法の一例を示す図である。

最適化装置１０は、図１の最適化装置１０ａと同様に、記憶部１１と処理部１２を有する。記憶部１１は前述の記憶部１１ａと、処理部１２は前述の処理部１２ａと、それぞれ同様のハードウェアにて実現できる。

記憶部１１は、前述の記憶部１１ａと同様の機能を有する。一方、処理部１２は、前述の処理部１２ａとは異なる以下に示すような機能を有する。
処理部１２は、レプリカが有する複数の状態変数のうちの第１状態変数の値を更新した場合における、前述の状態空間１内でのそのレプリカと、Ｍ個のレプリカからそのレプリカを除いたレプリカ群のうちの一部である他のレプリカとの距離の変化に応じた相互作用の強さの変化量を特定する。すなわち、第１の実施の形態に係る最適化方法では、全てのレプリカ間の相互作用が考慮されるのではなく、一部のレプリカ間の相互作用を考慮した処理が行われる。例えば、図２のように、レプリカ２，３間の相互作用とレプリカ３，４の相互作用は考慮されるが、レプリカ２，４間の相互作用は考慮されない。つまり、レプリカ２，４間には相互作用が与えられない。

各レプリカに対して相互作用を与えるレプリカの選択方法については、例えば、以下の４つがある。以下簡単に説明する。
１つ目の選択方法は、各レプリカに与えられている各レプリカを識別する識別情報であるレプリカ番号に基づいて、Ｍ個のレプリカに周期的に相互作用を与えるものである。この方法では、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカは、ｌ±ｓのレプリカ番号の範囲のものに限られる。つまり、ｌとの差がｓである範囲に含まれるレプリカ番号が与えられているレプリカが、レプリカ番号＝ｌのレプリカに対して相互作用を与えるものとなる。この方法では、相互作用の強さは、比較例と異なり、後述の式（１５）または式（１６）の代わりに、後述の式（２７）のように定義される。

２つ目の方法は、各レプリカに与えられたレプリカ番号に基づいてＭ個のレプリカを複数のグループにグループ分けし、異なるグループ間に属するレプリカ間についてだけ、相互作用を与えるものである。この方法では、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカは、レプリカ番号＝ｌのレプリカが属するグループとは異なる各グループの代表レプリカに限られる。この方法では、相互作用の強さは、比較例と異なり、後述の式（１５）または式（１６）の代わりに、後述の式（２８）のように定義される。

３つ目の方法は、相互作用を適用するレプリカの範囲を動的に決定する方法である。この方法では、１つ目の方法と同じように、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカが、ｌ±ｓ_ｔのレプリカ番号の範囲のものとなるが、このｓ_ｔは動的に変化する。状態遷移を繰り返す処理が行われるたびに、レプリカ番号＝ｌのレプリカとｌ±ｓ_ｔのレプリカ番号の範囲の各レプリカとの距離の平均値が計算され、その平均値と２つの閾値（Ｄ_１，Ｄ_２（Ｄ_１＜Ｄ_２））との比較結果に基づいて、ｓ_ｔが減少または増加する。例えば、斥力の相互作用を発生させる場合、上記距離の平均値が、Ｄ_１よりも小さければｓ_ｔは＋１され、Ｄ_２よりも大きければｓ_ｔは－１される。引力の相互作用を発生させる場合はこの逆となる。この方法では、相互作用の強さは、比較例と異なり、後述の式（１５）または式（１６）の代わりに、後述の式（２９）のように定義される。

４つ目の方法は、相互作用を適用するレプリカの範囲をランダムに決定する方法である。この方法では、状態遷移を繰り返す処理が行われるたびに、レプリカ番号＝ｌのレプリカに対し、他の各レプリカについて、所定の確率ｐで相互作用を与えるレプリカとして採用する。この方法では、相互作用の強さは、比較例と異なり、後述の式（１５）または式（１６）の代わりに、後述の式（３１）のように定義される。

これらの方法を用いた例については、後述の第２の実施の形態において説明する。
第１の実施の形態の最適化方法におけるその他の処理については、比較例の最適化方法の処理と同じである。すなわち、処理部１２は、前述のように、ΔＧ＝Ｇ（ｘ^ｌ［ｊ_０］）－Ｇ（ｘ^ｌ）と表される相互作用の強さの変化量を特定する。そして、処理部１２は、提案確率（ｇ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）と、受け入れ確率（ａ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）とに基づいて、確率的に第１状態変数（例えばｊ_０番目の状態変数）の値を更新した場合における、第１状態変数の値を更新するか否かを決定する。処理部１２は、第１状態変数の値を更新すると決定した場合、第１状態変数の値を更新した後のレプリカの複数の状態変数の値に基づいて、目的関数の値を計算する。また処理部１２は、記憶部１１内におけるレプリカの第１状態変数の値を更新する。そして処理部１２は、複数のレプリカそれぞれの複数の状態変数の一の状態変数の値の更新を繰り返し、目的関数の値が所定の条件を満たしたときの複数の状態変数の値を出力する。

このようにして、第１の実施の形態の最適化方法では、一部のレプリカ間に相互作用を与えたレプリカの状態遷移により、解の探索が行われる。すなわち最適化装置１０は、レプリカ間に相互作用を与えたことで、状態空間１内を網羅的に探索することができる。前述のように、全てのレプリカ間の相互作用を考慮しなくても、局所解からの脱出効果が得られる場合が多いためである。

このように、最適化装置１０は、全てのレプリカ間に相互作用を与えるのではなく、一部のレプリカ間に相互作用を与えるものであるため、比較例の最適化装置１０よりも計算量を抑えることができる。

〔第２の実施の形態〕
次に第２の実施の形態について説明する。第２の実施の形態は、目的関数の値が最小となる各状態変数の値の組み合わせを計算するイジングマシンを用いたシステムの例である。なお第２の実施の形態におけるイジングマシンは、第１の実施の形態に示した最適化装置１０の一例である。イジングマシンでは、求解対象の問題をイジングモデルで表し、そのイジングモデルのエネルギーが最小値となるビットの値の組み合わせを探索する。イジングモデルのエネルギーを計算する式（ハミルトニアン）が、目的関数である。

図３は、第２の実施の形態のシステム構成の一例を示す図である。サーバ１００には、ネットワーク２０を介して端末装置３１，３２，・・・が接続されている。端末装置３１，３２，・・・は、組み合わせ最適化問題の求解を依頼するユーザが使用するコンピュータである。サーバ１００は、端末装置３１，３２，・・・から組み合わせ最適化問題の求解の依頼を受け付け、組み合わせ最適化問題に対応するイジングモデルのエネルギー関数であるハミルトニアンを生成する。サーバ１００には、イジングマシン３００の制御装置２００が接続されている。サーバ１００は、生成したハミルトニアンを用いてエネルギーの最小値の探索要求を制御装置２００に入力する。

制御装置２００は、イジングマシン３００を制御し、サーバ１００から入力された探索要求に応じて、エネルギーの最小値の解探索を行う。例えば制御装置２００は、各ニューロンについての結合先のニューロンのｉｄを、結合先情報としてイジングマシン３００に送信する。また、制御装置２００は、ローカルフィールドの初期値（例えばバイアス係数）や、値が０ではない重み係数、アニーリング条件などについてもイジングマシン３００に送信する。

イジングマシン３００は、制御装置２００からの制御に基づいて、デジタル回路を用いたイジングモデルの状態遷移のシミュレーションを行い、エネルギーの最小値を探索する。

図４は、サーバのハードウェアの一例を示す図である。サーバ１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。

メモリ１０２は、サーバ１００の主記憶装置として使用される。メモリ１０２には、プロセッサ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、メモリ１０２には、プロセッサ１０１による処理に利用する各種データが格納される。メモリ１０２としては、例えばＲＡＭ（Random Access Memory）などの揮発性の半導体記憶装置が使用される。

バス１０９に接続されている周辺機器としては、ストレージ装置１０３、グラフィック処理装置１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７およびネットワークインタフェース１０８がある。

ストレージ装置１０３は、内蔵した記録媒体に対して、電気的または磁気的にデータの書き込みおよび読み出しを行う。ストレージ装置１０３は、コンピュータの補助記憶装置として使用される。ストレージ装置１０３には、ＯＳのプログラム、アプリケーションプログラム、および各種データが格納される。なお、ストレージ装置１０３としては、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）を使用することができる。

グラフィック処理装置１０４には、モニタ２１が接続されている。グラフィック処理装置１０４は、プロセッサ１０１からの命令に従って、画像をモニタ２１の画面に表示させる。モニタ２１としては、有機ＥＬ（Electro Luminescence）を用いた表示装置や液晶表示装置などがある。

入力インタフェース１０５には、キーボード２２とマウス２３とが接続されている。入力インタフェース１０５は、キーボード２２やマウス２３から送られてくる信号をプロセッサ１０１に送信する。なお、マウス２３は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。

光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク２４に記録されたデータの読み取り、または光ディスク２４へのデータの書き込みを行う。光ディスク２４は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク２４には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ－ＲＡＭ、ＣＤ－ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ－Ｒ（Recordable）／ＲＷ（ReWritable）などがある。

機器接続インタフェース１０７は、サーバ１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置２５やメモリリーダライタ２６を接続することができる。メモリ装置２５は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ２６は、メモリカード２７へのデータの書き込み、またはメモリカード２７からのデータの読み出しを行う装置である。メモリカード２７は、カード型の記録媒体である。

ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。ネットワークインタフェース１０８は、例えばスイッチやルータなどの有線通信装置にケーブルで接続される有線通信インタフェースである。またネットワークインタフェース１０８は、基地局やアクセスポイントなどの無線通信装置に電波によって通信接続される無線通信インタフェースであってもよい。

サーバ１００は、以上のようなハードウェアによって、第２の実施の形態の処理機能を実現することができる。なお、制御装置２００も、サーバ１００と同様のハードウェアにより実現することができる。

図５は、イジングマシンの一例を示す図である。イジングマシン３００は、ニューロン回路３１１，３１２，…，３１ｎ、制御回路３２０、およびメモリ３３０を有する。
ニューロン回路３１１～３１ｎはそれぞれ、自身以外の複数の他のニューロン回路との接続の有無を示す複数の重み係数の値と、複数の他のニューロン回路の複数の出力信号との積の総和に基づく第１の値を算出する。そしてニューロン回路３１１～３１ｎそれぞれは、第１の値にノイズ値を加算した第２の値と閾値との比較結果に基づき、０または１のビット値を出力する。複数のレプリカを用いた解探索を行う場合、複数のニューロン回路を用いて１つのレプリカの解探索が行われる。

制御回路３２０は、制御装置２００から供給される情報に基づいて、イジングマシン３００の初期設定処理などを行う。また、制御回路３２０は、レプリカ交換を行う場合、２つのレプリカ間の温度パラメータの値の交換の有無を判断し、温度パラメータの値を交換する場合、各レプリカの解探索を行うニューロン回路に入力する温度パラメータの値を更新する。

さらに、制御回路３２０は、更新対象ニューロンを決定する処理が所定回数だけ繰り返された後、メモリ３３０に保持されている１レプリカの状態変数に対応する各ニューロンのビット値を取得し、最適化問題に対する解として制御装置２００に送信する。

制御回路３２０は、例えばＡＳＩＣやＦＰＧＡなどの特定用途の電子回路にて実現できる。なお、制御回路３２０は、ＣＰＵやＤＳＰなどのプロセッサであってもよい。その場合、プロセッサは、図示しないメモリに記憶されたプログラムを実行することで、上記の処理を行う。

メモリ３３０は、例えば各ニューロンのビット値を保持する。メモリ３３０は、例えばレジスタやＲＡＭなどによって実現できる。メモリ３３０には、エネルギーの最小値や最小値が得られたときの各ニューロンのビット値を保持しておくこともできる。この場合、制御回路３２０は、更新対象ニューロンを決定する処理が所定回数、繰り返されたのちに、エネルギーの最小値や最小値が得られたときの各ニューロンのビット値をメモリ３３０から取得して、制御装置２００に送信してもよい。

なお第１の実施の形態に示した最適化装置１０も、図４に示したイジングマシン３００と同様のハードウェアにより実現することができる。
次に、求解の対象となるイジング型の最小値求解問題（イジング型問題）について説明する。イジング型問題はイジングモデルで表される。

図６は、イジングモデルの模式図である。イジングモデル３０は、格子状に複数のビット３１が配置される。ビット３１は、磁石を模したものであり、スピンとも呼ばれる。隣接するビット間には相互作用が働く。相互作用の大きさは、重み係数で表される。イジングモデル３０のエネルギーは、次の式（１）で表される。

右辺の１項目は、Ｎ個の状態変数の全組み合わせについて、漏れと重複なく、２つの状態変数の値（０または１）と重み係数との積を積算したものである。ｘ_ｉはｉ番目の状態変数、ｘ_ｊはｊ番目の状態変数を表し、Ｗ_ｉｊは、ｘ_ｉとｘ_ｊとの結合の強さを示す重み係数である。右辺の２項目は、各状態変数のそれぞれについてのバイアス係数（ｂ_ｉ）とｘ_ｉとの積の総和を求めたものである。Ｗ_ｉｊが正の場合、ｘ_ｉとｘ_ｊとが同じ値となるような相互作用が働く。またＷ_ｉｊが負の場合、ｘ_ｉとｘ_ｊとが異なる値となるような相互作用が働く。なおＷ_ｉｊ＝Ｗ_ｊｉであり、Ｗ_ｉｉ＝０である。

最小値求解問題は、式（１）で与えられるエネルギーの最小値を求める問題である。イジングマシン３００は、このような最小値求解問題についてＭＣＭＣを用いて解く。例えばイジングマシン３００は、ビットを１つ反転した場合のエネルギー変化を計算する。ｉ番目のビットを反転させた場合「ｘ_ｉ→ｘ_ｉ´（δｘ_ｉ＝ｘ_ｉ´－ｘ_ｉ）」、エネルギーの変化値は式（２）で表される。

式（２）の右辺の括弧内の式は、ｉ番目のビットの局所場（総入力）を表している。出力変化分δｘ_ｉと局所場の符号が一致すればエネルギーは減少する。イジングマシン３００は、エネルギーの変化値ΔＥ_ｉの増減に応じて、ｉ番目のビットの反転を受け入れるか否かを決定する。なお、式（２）は、１ビットのみを反転させる場合にのみ正しい式である。

エネルギー増分を示す式（２）は、以下のように書き換えることができる。

ｈ_ｉはｉ番目のビットの局所場である。ｊ番目のビットｘ_jが反転したときのｉ番目のビットの局所場ｈ_ｉの変化分δｈ_ｉ ^（ｊ）は、以下の式（５）で表される。

局所場ｈ_ｉを保存するレジスタを用意し、ｊ番目のビットが反転したときに、式（５）に示す値を、保存しておいた局所場ｈ_ｉに加算することで、常に正しいｈ_ｉが得られる。
以上のような計算により、ｉ番目のビットの反転させた場合のエネルギーの増分を求めることができる。イジングマシン３００は、得られたエネルギーの増分に基づいて、ｉ番目のビットの反転を受け入れるか否かを判断する。例えばイジングマシン３００は、メトロポリス法に従って、ビットの反転を受け入れるか否かを判断する。メトロポリス法に従う場合、エネルギーの増分が負（エネルギーが減少する）であれば、ビットの反転が受け入れられる。エネルギーの増分が正（エネルギーが増加する）であれば、エネルギーの増分に応じた確率によって、ビットの反転が受け入れられるか否かが判断される。

エネルギーの増分が正の場合にビットの反転が受け入れられる確率は、温度パラメータを用いて調整することができる。例えばイジングマシン３００は、温度パラメータの値が大きいほど、エネルギーの増分が正の場合にビットの反転を受け入れる確率を高くする。これにより、温度パラメータの値を大きくすることで、イジングモデルのエネルギーの状態が局所解から抜け出す可能性を高くすることができる。

温度パラメータをＴとしたとき、逆温度β＝１／Ｔとする。例えばイジングマシン３００は、エネルギーの変化値ΔＥ_ｉｊと逆温度βを用いて、ｉ番目の状態変数の状態遷移の受け入れ確率を以下の式（６）により決定することで、確率的探索を行う。

式（６）における関数ｆ（ｘ）は、メトロポリス法では以下の式（７）となる。

なお温度パラメータの値が大きいと、局所的な探索が難しくなる。そこでイジングマシン３００は、例えば温度パラメータの値が異なる複数のレプリカを用いて解探索を行う。この場合、イジングマシン３００はレプリカ交換を行ってもよい。

図７は、レプリカ交換の一例を示す図である。レプリカ交換では、複数のレプリカが用いられる。レプリカは、求解対象の問題の状態変数の集合のコピーである。イジングマシン３００は、それぞれのレプリカの温度パラメータに異なる値を設定する。図７の例では、４つのレプリカそれぞれに、Ｔ_１，Ｔ_２，Ｔ_３，Ｔ_４の温度パラメータが設定されている（Ｔ_１＜Ｔ_２＜Ｔ_３＜Ｔ_４）。

イジングマシン３００は、複数のレプリカそれぞれについて、ＭＣＭＣにより状態変化させる。そしてイジングマシン３００は、温度パラメータの値で並べたときに隣り合うレプリカ間で、所定の確率に従って温度パラメータの値を交換する。すると、各レプリカは温度軸方向にランダムウォークする。レプリカがランダムウォークをすることで、局所解に嵌まっても高温側に移動したときに局所解から脱出できる可能性がある。またレプリカが低温側に移動すると、局所的な探索を行うこともできる。

レプリカ交換のように多くのレプリカを用いて集団探索を行うことで、モンテカルロ法による解探索を高速化することができる。しかし単に複数のレプリカを用いて集団探索をしただけでは、複数のレプリカが同じ局所解に留まってしまい、十分に広く状態空間を探索できないという問題を解決できない。例えばイジングモデルの状態変数（ビット）の数がＮのとき、状態空間には２^Ｎ個の状態が存在する。そのため状態変数の数が多くなると、実用上可能な数のレプリカで探査しても、集団探索の利益を享受するのは困難である。

そこでイジングマシン３００では、一部のレプリカ間の距離に応じた相互作用を利用して、レプリカの状態を遷移させることで、状態空間内の効率的な探索を行う。これにより、複数のレプリカを用いた集団探索による解探索性能が向上する。

例えばレプリカ交換を行えば状態空間内を広範囲に探索できるが、レプリカ間の相互作用を考慮しない場合、各レプリカは、そのときの温度パラメータの値に応じて、独立にビットフリップ（マルコフ連鎖）を行うだけである。レプリカ間の相互作用を利用すれば、個々のレプリカのマルコフ連鎖において、複数のレプリカが同じ局所解に同時に留まることを抑止可能である。

また１－ｂｉｔフリップの場合、遷移先の候補の選び方としてＮ個のビットを等確率で選択すると、遷移確率はエネルギーの変化値ΔＥ_ｉのみで決定される。この場合、どのレプリカもエネルギー勾配のみに従って状態変化するため、同じ道筋を辿ってしまい、十分に広く状態空間を探索できていない可能性が大きい。さらに、どのレプリカも同じ局所解（すべてのビットｉでΔＥ_ｉ＞０）に嵌まってしまったときに脱出するのも困難である。

イジングマシン３００では、遷移を受け入れるか否かの受容確率の計算にメトロポリス法ではなく、メトロポリス・ヘイスティングス法を用いる。これにより、レプリカ間の相互作用による影響を適切な形で計算に組み込むことができる。

例えば現在の状態ｘから次の状態Ｘ′を提案する確率をｇ（Ｘ→Ｘ′）とし、この状態遷移が受け入れられる確率をＡ（Ｘ→Ｘ′）とする。状態Ｘから状態Ｘ′に遷移する確率Ｗ（Ｘ→Ｘ′）は、以下の式（８）で得られる。

目的の確率分布（例えばギブス分布）を表す関数（目的関数）をπ（Ｘ）とすると詳細つり合いの条件は以下の通りとなる。

式（１０）から、詳細つり合いを満たす受け入れ確率は、式（１１）の通りとなる。

メトロポリス・ヘイスティングス法を適用した場合、受け入れ確率は以下の式（１２）で与えられる。

この受け入れ確率は、提案確率が非対称でｇ（Ｘ→Ｘ′）≠ｇ（Ｘ′→Ｘ）の場合でも詳細つり合いの条件は満たされる。また提案確率が対称でｇ（Ｘ→Ｘ′）＝ｇ（Ｘ′→Ｘ）の場合は式（１３）のようなメトロポリスの受け入れ確率が得られる。

ここで１－ｂｉｔフリップを考えるとき、レプリカ間の相互作用を考慮しない場合、Ｎ個のビットが等確率で反転の候補として選択され、提案確率は式（１４）となる。

なお、メトロポリス・ヘイスティングス法は提案確率で示される提案分布が非対称な場合に対応している。そのため、提案分布の決め方に自由度がある。そこでイジングマシン３００は、提案確率にレプリカ間相互作用を導入する。

例えばイジングマシン３００は、離散空間である状態空間に対して適切な距離を定義し、レプリカ間の距離を定める。イジングマシン３００は、レプリカ間距離を用いてレプリカ相互作用を決め、メトロポリス・ヘイスティングス法における遷移候補先の分布（提案分布）を定義し、受け入れ確率の計算に組み込む。

レプリカ間距離の一例として、２つのレプリカの状態のハミング距離（或いはその単調増加関数）がある。イジングマシン３００は、すべてのレプリカ間のハミング距離を定義し、それによりレプリカ相互作用を導入する。

式（１４）に示すような提案分布は１／Ｎの等確率で遷移候補が選択されていたため、規格化定数はＮ（それぞれの遷移先の重みが共通で１ということ）である。レプリカ間の相互作用を導入した場合は遷移候補の重みが異なり、規格化定数が遷移前の現在状態に依存する。イジングマシン３００は、規格化定数も計算しなければならないが、相互作用がハミング距離の一次式の場合は、差分計算（累積計算）によって容易に計算できる。

以下、レプリカ間距離を考慮した提案確率の計算方法について具体的に説明する。まず提案確率の一般系を以下のように定義する。
Ｍ個（Ｍは１以上の整数）のレプリカからなる系を考える。１番目のレプリカの状態変数をｘ^ｌ＝（ｘ_１ ^ｌ，ｘ_２ ^ｌ，・・・，ｘ_Ｎ ^ｌ），ｘ_ｊ ^ｌ∈｛０，１｝とする。２つのレプリカｘ^ｌとｘ^ｋの距離（の増加関数）をｄ（ｘ^ｌ，ｘ^ｋ）とし、相互作用のエネルギーをＧ（ｘ）と与える。相互作用のエネルギーは、全レプリカ間の相互作用を考慮した場合、例えば式（１５）または式（１６）のようにいくつか定義することができる。

ここでγは実数の定数である。γが正の値であれば引力的な相互作用、γが負の値であれば斥力的な相互作用とみなせる。
式（１５）または式（１６）に含まれるレプリカ間の距離の計算回数は、各レプリカについてＭ回であるため、イジングマシン３００全体では、Ｍ^２回となり、レプリカ数が増えるほど計算量が大幅に増大する。このため、イジングマシン３００では、Ｍ個のレプリカのうちの一部のレプリカ間に相互作用を与えて処理を行う。各レプリカに対して相互作用を与えるレプリカの選択方法については、例えば、後述の４つの方法があり、４つの方法のそれぞれにおいて定義されるＧ（ｘ）が異なる。各方法で用いられるＧ（ｘ）については後述する。

Ｇ（ｘ）を用いて提案確率をｇ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）と与える。ｘ^ｌ［ｊ_０］はｘ^ｌのｊ_０番目のビットをフリップした状態を表す。提案確率は、具体的には式（１７）または式（１８）のように定義することができる。

ここでＺ（ｘ^ｌ）は規格化定数であり、計算方法は後述する。
レプリカ間距離としてハミング距離の１次関数を用いた場合、レプリカ間距離は式（１９）で定義できる。

この場合、ΔＧ＝Ｇ（ｘ^ｌ［ｊ_０］）－Ｇ（ｘ^ｌ）とｇ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）は以下のように計算できる。

このようにして、レプリカ間の相互作用を反映させた提案確率を計算することができる。次に、受け入れ確率の定義について説明する。
一般系の受け入れ確率ａ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）は、メトロポリス基準を採用することにすると、以下のように定義できる。

すると遷移確率はＷ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）＝ｇ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）×ａ（ｘ^ｉ→ｘ^ｉ［ｊ_０］）となる。従って、これらの計算に用いられる量はΔＥ、ΔＧ、規格化定数Ｚの３つである。

ここで、レプリカ間距離としてハミング距離の一次関数を用いた場合を例にとり、規格化定数Ｚの計算方法について説明する。規格化定数Ｚ（ｘ^l）は、提案候補がレプリカｌにおいて１つのビットをフリップしただけの状態なので、その総和として以下の式（２３）で計算される。

式（２３）のまま規格化定数を計算しようとするとすべてのスピンの数だけ指数関数の和を計算することとなり、計算量が膨大となる。そこでイジングマシン３００は、１－ｂｉｔフリップであることに基づいて、差分計算（累積計算）を行うことで、計算量を抑止する。レプリカｌのｊ番目のビットだけをフリップした場合の規格化定数の差分は以下の式（２４）のようになる。

イジングマシン３００は、式（２４）の右辺を計算することで求めた規格化定数の差分を、ビットフリップ前の規格化定数に加算することで、ビットフリップ後の規格化定数を求めることができる。なおイジングマシン３００は、ビットフリップが受け入れられた場合、そのときの規格化定数をレジスタまたはメモリに保存し、次回のビットフリップにおける規格化定数の算出に利用する。

次に、レプリカ間距離としてハミング距離の一次関数を用いた場合を例にとり、ΔＧの計算方法について説明する。ΔＧの計算は、一般にはレプリカ間の距離（あるいは距離の増加関数）の差分計算になる。単純な差分計算では、遷移前後でのレプリカ間のハミング距離を記憶しておかなければならない。しかし、具体的に距離（あるいは距離の増加関数）の形が分かれば、差分計算を行うことで、式（２５）、式（２６）に示すように、現在の状態だけに依存した量に書き換えることができる。

式（２６）において、ｘ_ｊ０（ｌ）（ｘはチルダ付き）は、ｘ_ｊ０（ｌ）（ｘはチルダ付き）＝（ｘ_ｊ ^ｌ，ｘ_ｊ ^ｌ，・・・，ｘ_ｊ ^ｌ）というビット列である。またｘ_ｊ０（ｘはチルダ付き）は、ｘ_ｊ０（ｘはチルダ付き）＝（ｘ_ｊ ^１，ｘ_ｊ ^２，・・・，ｘ_ｊ ^Ｍ）というビット列のベクトルである。

式（２６）を用いれば、レプリカ間距離をハミング距離の一次関数とした場合には、新たに導入したビット列のベクトル間のハミング距離だけで、ΔＧを記述することができる。従って、そのハミング距離だけを更新すればよい。

なお、これまでの説明では、１－ｂｉｔフリップの場合を想定しているが、１回の状態遷移で複数のビットをフリップする場合もある。例えばＯｎｅ－Ｈｏｔ制約がある問題を求解する場合である。

Ｏｎｅ－Ｈｏｔ制約は「ある変数の組の中で、値が１の変数は１個だけ」という制約である。この制約は、二次割り当て問題（ＱＡＰ：Quadratic Assignment Problem）、運搬経路問題（ＶＲＰ：Vehicle Routing Problem）などの様々な問題に適用される。

図８は、Ｏｎｅ－Ｈｏｔ制約下での１－ｂｉｔフリップの例を示す図である。図８の例では、イジングモデルの状態変数を示すビットが、４ビットずつのグループに分けられている。Ｏｎｅ－Ｈｏｔ制約では、同一グループに属するビットのうち１ビットのみが「１」であることを許される。このようなＯｎｅ－Ｈｏｔ制約の元で１－ｂｉｔフリップを行うと、１回の状態遷移では１ビットのみが反転し、制約違反の状態となる。１－ｂｉｔフリップをもう１回行うと、制約条件を満足することができる。

このようにＯｎｅ－Ｈｏｔ制約を有する問題を求解する場合、１－ｂｉｔフリップでは効率が悪い。そこで、イジングマシン３００は、１回の状態遷移で複数のビットをフリップさせることができる。

Ｏｎｅ－Ｈｏｔ制約には、１－Ｗａｙ １－Ｈｏｔ（１Ｗ１Ｈ）と２－Ｗａｙ １－Ｈｏｔ（２Ｗ１Ｈ）とがある。１Ｗ１Ｈは、１つの手段でビットをグルーピングしたときの各グループ内で、値が「１」のビットは１つのみという制約である。図８に示した例は１Ｗ１Ｈであり、１回の状態遷移で２つのビットをフリップさせることで、制約を満たしたままの状態遷移が可能となる。

２Ｗ１Ｈでは、２つの手段でビットがグルーピングされる。この場合、各ビットは、生成手段が異なる２つのグループに属する。そして２Ｗ１Ｈでも、各グループ内で値が「１」のビットは１つのみという制約がある。

図９は、２Ｗ１Ｈの制約を説明する図である。図９には、Ｎ個のビットをｎ×ｎ（ｎは１以上の整数）の正方行列の要素としている。Ｎ＝ｎ²である。２Ｗ１Ｈでは、各行および各列のビットの値の和は「１」であるという制約がある。すなわち同じ行のビットのうちの１つのビットだけ値が「１」であり、同じ列のビットのうちの１つのビットだけ値が１である場合に制約が満たされる。２Ｗ１Ｈの制約がある場合、１回の状態遷移で４つのビットをフリップさせることで、制約を満たしたままの状態遷移が可能となる。

ここでｍ＝１，２，・・・，Ｎとしたとき、１－ｂｉｔフリップ、１Ｗ１Ｈにおける２－ｂｉｔフリップ、２Ｗ１Ｈにおける４－ｂｉｔフリップそれぞれの、状態遷移、エネルギーの変化値ΔＥ、局所場の更新量Δｈは、以下のように表される。

＜１－ｂｉｔフリップ＞
・状態遷移 ：ｘ_ｉ→ｘ_ｉ＋Δｘ_ｉ
・エネルギーの変化値：ΔＥ_ｉ＝－ｈ_ｉ・Δｘ_ｉ
・局所場更新量 ：Δｈ_ｍ＝Ｗ_ｍｉ・Δｘ_ｉ
＜１Ｗ１Ｈ（２－ｂｉｔフリップ）＞
・状態遷移：ｘ_ｉ：１→０，ｘ_ｊ：０→１
・エネルギーの変化値：ΔＥ_ｊ＝ｈ_ｉ－ｈ_ｊ
・局所場更新量：Δｈ_ｍ＝－Ｗ_ｍｉ＋Ｗ_ｍｊ
＜２Ｗ１Ｈ（４－ｂｉｔフリップ）＞
・状態遷移：ｘ_ｉ：１→０，ｘ_ｊ：０→１，ｘ_ｋ：０→１，ｘ_ｌ：１→０
・エネルギーの変化値：ΔＥ_ｊ＝（ｈ_ｉ＋ｈ_ｌ）－（ｈ_ｊ＋ｈ_ｋ）－（Ｗ_ｉｌ＋Ｗ_ｊｋ）
・局所場更新量：Δｈ_ｍ＝Ｗ_ｍｊ＋Ｗ_ｍｋ－（Ｗ_ｍｉ＋Ｗ_ｍｌ）
いずれの制約を適用するのかは、例えば、ユーザが問題の求解を指示する際に、ユーザによって指定される。イジングマシン３００は、指定された制約に応じたΔＥを計算し、レプリカ間の距離に応じた遷移確率で１または複数のビットを反転させる。

次に、レプリカ間距離を考慮したイジングマシン３００による解探索機能について説明する。
図１０は、イジングマシンの解探索機能の一例を示す図である。イジングマシン３００は、データ受け取り部３４０、解探索エンジン３５０、および解出力部３６０を有している。データ受け取り部３４０と解出力部３６０とは、図５に示した制御回路３２０によって実現される機能である。解探索エンジン３５０は、図５に示した制御回路３２０が、ニューロン回路３１１，３１２，・・・，３１ｎとメモリ３３０とを制御することで実現する機能である。

データ受け取り部３４０は、制御装置２００から探索対象の問題の求解に用いる情報を受け取る。例えばデータ受け取り部３４０は、温度、レプリカ数、レプリカ相互作用の大きさ、iteration数（状態遷移の反復回数）、初期状態などのパラメータを取得する。またデータ受け取り部３４０は、求解対象の問題を表すイジングモデルの重み係数を要素とする重み行列（２次式の係数）、バイアス行列（１次式の係数）、定数項、１－Ｈｏｔ制約のグループ情報などのデータを取得する。さらに、データ受け取り部３４０から、レプリカ間の相互作用を適用する範囲を決めるための後述のパラメータを取得する。データ受け取り部３４０は、受け取った情報を解探索エンジン３５０に送信する。

解探索エンジン３５０は、複数のレプリカを用いて、エネルギーが最小となる解を探索する。そのために解探索エンジン３５０は、レプリカ保存部３５１と複数のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎとを有する。レプリカ保存部３５１は、例えば図５に示したメモリ３３０を利用して実現される。複数のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、イジングモデルに含まれるビットごとのニューロン回路を利用して実現される。

レプリカ保存部３５１は、レプリカの状態を記憶する。例えばレプリカが順番に更新されていくが、レプリカ間相互作用の計算には、更新前のレプリカの状態が使用される。そこでレプリカ保存部３５１が、更新前のレプリカの状態を記憶する。レプリカの状態は、状態変数に対応するビットの値、および温度パラメータなどのパラメータの値で表される。

各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎは、それぞれがレプリカによる解探索を行う。例えば各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎは、レプリカ保存部３５１を介して互いのレプリカの状態を示す情報をやり取りしながらレプリカ間相互作用を計算し、解の探索を行う。

図１１は、解探索エンジンにおける処理の一例を示す図である。例えばレプリカ解探索部３５２ａは、重み係数（Ｗ_ｉｊ）を記憶している。レプリカ解探索部３５２ａは、重み係数（Ｗ_ｉｊ）と現在の各ビットの値（ｘ_１ ^１,ｘ_２ ^１，・・・，ｘ_Ｎ ^１）とを用いて、式（４）に基づいて局所場（ｈ_１，ｈ_２，・・・，ｈ_Ｎ）を計算する。次にレプリカ解探索部３５２ａは、式（２６）に基づいて、各ビットがフリップした場合のレプリカ間の相互作用のエネルギーの差分（ΔＧ_１，ΔＧ_２，・・・，ΔＧ_Ｎ）を計算する。この際、レプリカ解探索部３５２ａは、レプリカ保存部３５１から、他のレプリカ（Ｍ個のレプリカから自身が解探索を担当するレプリカを除いたレプリカ群のうち一部）の状態を示す情報（各ビットの値）を取得し、他のレプリカとの距離を計算し、その計算結果を用いてレプリカ間の相互作用のエネルギーの差分を計算する。

さらにレプリカ解探索部３５２ａは、局所場（ｈ_１，ｈ_２，・・・，ｈ_Ｎ）の値を用いて、エネルギーの変化値（Ｅ_１，Ｅ_２，・・・，Ｅ_Ｎ）を計算する。なおエネルギーの変化値の計算式は、１－ｂｉｔフリップなのか１Ｗ１Ｈなのか２Ｗ１Ｈなのかによって異なる。例えば１－ｂｉｔフリップであれば、エネルギーの変化値は「ΔＥ_ｉ＝－ｈ_ｉ・Δｘ_ｉ」である。１Ｗ１Ｈ（２－ｂｉｔフリップ）であれば、エネルギーの変化値は「ΔＥ_ｊ＝ｈ_ｉ－ｈ_ｊ」である。２Ｗ１Ｈ（４－ｂｉｔフリップ）であれば、エネルギーの変化値はΔＥ_ｊ＝（ｈ_ｉ＋ｈ_ｌ）－（ｈ_ｊ＋ｈ_ｋ）－（Ｗ_ｉｌ＋Ｗ_ｊｋ）である。

レプリカ解探索部３５２ａは、エネルギーの変化値ΔＥから正のオフセット値Ｅ_ｏｆｆを減算する。オフセット値Ｅ_ｏｆｆには、フリップするビットが選択できなかった場合に、所定の値が加算される。オフセット値Ｅ_ｏｆｆの増加は、フリップするビットが選択されるまで繰り返される。このように、オフセット値Ｅ_ｏｆｆが増加することで、レプリカのエネルギーが極小値に留まる時間が短縮される。なお、オフセット値Ｅ_ｏｆｆの初期値は、例えば「０」とする。

レプリカ解探索部３５２ａは、各ビットをフリップさせた場合のエネルギーの変化値ΔＥ（オフセット値Ｅ_ｏｆｆが「０」以外の場合にはオフセット値Ｅ_ｏｆｆを減算後の値）に基づいて、フリップするビット（更新ビット）を選択する。更新ビットの選択方法には、様々な方法がある（図２５～図２８参照）。更新ビットの選択方法によっては、更新ビットの選択において、いずれのビットの更新の受け入れも棄却され、更新ビットが選択できないことがあり得る。レプリカ解探索部３５２ａは、例えば更新ビットが選択できなかった場合、オフセット値Ｅ_ｏｆｆの値を増加させ、再度、更新ビットの選択を行う。

レプリカ解探索部３５２ａは、更新ビットが選択できた場合、更新ビットの値をフリップし、更新後のレプリカの状態「ｘ_１ ^１，ｘ_２ ^１，・・・，ｘ_Ｎ ^１」を生成する。
レプリカ解探索部３５２ａ以外のレプリカ解探索部３５２ｂ，・・・，２５２ｎも、レプリカ解探索部３５２ａと同様に、更新後のレプリカの状態を生成する。

各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎが生成したレプリカの状態「ｘ_１ ^１，ｘ_２ ^１，・・・，ｘ_Ｎ ^１」、「ｘ_１ ^２，ｘ_２ ^２，・・・，ｘ_Ｎ ^２」、・・・、「ｘ_１ ^Ｎ，ｘ_２ ^Ｎ，・・・，ｘ_Ｎ ^Ｎ」は、レプリカ保存部３５１で保持される。各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎは、レプリカ保存部３５１を参照することで、次回の状態更新時に、レプリカ間の相互作用のエネルギーの差分を算出することができる。

以下、解探索エンジン３５０による解探索の手順について詳細に説明する。
解探索の手順は、各レプリカに対して相互作用を与えるレプリカの選択方法によって異なる。

（相互作用を与えるレプリカの１つ目の選択方法）
１つ目の選択方法は、各レプリカに与えられたラベル（レプリカ番号）に基づいて、Ｍ個のレプリカに周期的に相互作用を与えるものであり、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカは、ｌ±ｓのレプリカ番号の範囲のものに限られる。

１つ目の選択方法では、レプリカ番号＝ｌのレプリカに対して与えられる相互作用の強さは、式（１５）または式（１６）の代わりに、例えば、以下の式（２７）で定義できる。

すなわち、レプリカ番号＝ｌのレプリカと、ｌ±ｓのレプリカ番号の範囲のレプリカとの距離に基づいて、相互作用の強さが定義される。
図１２は、相互作用を与えるレプリカの１つ目の選択方法の一例を示す図である。

図１２では、１２個のレプリカに対して、１～１２のレプリカ番号が与えられており、このうち、各レプリカに相互作用を与えるレプリカの範囲（相互作用の適用範囲）を決めるｓが、ｓ＝２である例が示されている。図１２のように、レプリカ番号＝１のレプリカに対して相互作用を与えるレプリカは、１±２のレプリカ番号の範囲のレプリカとなる。この場合、１－２は負のレプリカ番号となってしまう。負のレプリカ番号が発生することを避けるために、図１２の例のように、ｌ＝１の１つ前はｌ＝１２（ｌ＝１２の次はｌ＝１）というように、レプリカ番号が循環していると考える。つまり、ｋ＝ｌ±ｓ ｍｏｄ Ｍが、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカの範囲であるものとする。すなわち、ｌ－ｓをＭで割った余りと、ｌ＋ｓをＭで割った余りの範囲に含まれるレプリカ番号が与えられているレプリカがレプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカの範囲となる。これにより、式（２７）において、ｋ＝Ｍ＋１の場合、ｋ＝１、ｋ＝－２の場合は、ｋ＝Ｍ－１となる。

図１３は、相互作用を与えるレプリカの１つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。以下、図１３に示す処理をステップ番号に沿って説明する。

［ステップＳ１００］解探索エンジン３５０は、レプリカ間の相互作用を適用する範囲を、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。相互作用の適用範囲は、前述のパラメータ（ｓ）によって決定される。ｓは、例えば、制御装置２００からデータ受け取り部３４０を介して解探索エンジン３５０に供給される。

［ステップＳ１０１］解探索エンジン３５０は、複数のレプリカの初期状態（各ビットの値、温度パラメータの値など）を、そのレプリカの割り当て先のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカの初期状態に基づいて、初期エネルギー、初期のレプリカ間距離、初期規格化定数などを計算する。

［ステップＳ１０２］解探索エンジン３５０は、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎにレプリカごとの解探索を実行させる。レプリカごとの解探索処理の詳細は後述する（図１４参照）。

［ステップＳ１０３］解探索エンジン３５０は、解探索の終了条件を満たしたか否かを判断する。例えば解探索エンジン３５０は、ステップＳ１０２の処理の繰り返し回数が所定回数に達した場合に、終了条件を満たすと判断する。解探索エンジン３５０は、終了条件を満たした場合、処理をステップＳ１０８に進める。また解探索エンジン３５０は、終了条件が満たされていない場合、処理をステップＳ１０４に進める。

［ステップＳ１０４］解探索エンジン３５０は、複数のレプリカを温度パラメータの値で並べたときに隣接するレプリカの組を選択する。
［ステップＳ１０５］解探索エンジン３５０は、選択したレプリカの組の温度交換の実施の有無を決定する。例えば解探索エンジン３５０は、レプリカ間のエネルギーの差と各レプリカの温度パラメータの値とに基づいて、メトロポリスヘイスティング基準により交換確率を求める。そして解探索エンジン３５０は、交換確率が１であれば温度交換を実施すると判断する。また解探索エンジン３５０は、交換確率が１未満であれば、例えば０から１までの間の乱数を生成し、乱数の値が交換確率以下であれば、温度交換を実施すると判断する。

［ステップＳ１０６］解探索エンジン３５０は、温度交換を実施すると決定した場合、選択したレプリカの組それぞれの温度パラメータの値を交換する。
［ステップＳ１０７］解探索エンジン３５０は、隣接するレプリカのすべての組を選択したか否かを判断する。解探索エンジン３５０は、未選択の組がある場合、処理をステップＳ１０４に進める。また解探索エンジン３５０は、すべての組が選択済みの場合、処理をステップＳ１０２に進める。

［ステップＳ１０８］解探索エンジン３５０は、エネルギーが最小となるレプリカの状態を、解として出力する。
このようにして、レプリカ交換を行いながら、複数のレプリカを用いた効率的な解探索が行われる。

次にレプリカごとの解探索処理について詳細に説明する。
図１４は、レプリカごとの解探索処理の手順の一例を示すフローチャートである。以下、図１４に示す処理をステップ番号に沿って説明する。

［ステップＳ１１０］解探索エンジン３５０内のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカについて、レプリカ間の相互作用のエネルギーの差分（ΔＧ_１，ΔＧ_２，・・・，ΔＧ_Ｎ）を計算する。レプリカ間の相互作用のエネルギーの差分の計算処理の詳細は後述する（図１５参照）。

［ステップＳ１１１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカについて、エネルギーの変化値（ΔＥ_１，ΔＥ_２，・・・，ΔＥ_Ｎ）を計算する。

［ステップＳ１１２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、反復回数をインクリメントする。
［ステップＳ１１３］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、所定回数だけ反復したか否かを判断する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、所定回数だけ反復した場合、レプリカごとの解探索処理を終了する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、反復回数が所定回数に達していなければ、処理をステップＳ１１４に進める。

［ステップＳ１１４］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビット選択処理を行う。更新ビット選択処理の詳細は後述する（図２５～図２８参照）。

［ステップＳ１１５］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビットが選択されたか否かを判断する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビットが選択されていない場合、処理をステップＳ１１４に進める。またレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビットが選択された場合、処理をステップＳ１１６に進める。

［ステップＳ１１６］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、レプリカに関する情報を更新する。例えばレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、選択されたビットの状態をフリップさせ、各ビットの局所場ｈ、レプリカのエネルギーＥ、他のレプリカとのレプリカ間距離ｄ、規格化定数Ｚを更新する。その後、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは処理をステップＳ１１０に進める。

次に、レプリカ間の相互作用のエネルギーの差分（ΔＧ_１，ΔＧ_２，・・・，ΔＧ_Ｎ）の計算処理について詳細に説明する。
図１５は、レプリカ間の相互作用のエネルギーの差分の計算手順の一例を示すフローチャートである。以下、図１５に示す処理をステップ番号に沿って説明する。

［ステップＳ１２０］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、自身が解探索を担当するレプリカと、他の何れかのレプリカとの間のハミング距離を計算する。相互作用を与えるレプリカの１つ目の選択方法では、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、解探索を担当するレプリカと、そのレプリカのレプリカ番号に対して±ｓの範囲のレプリカ番号のレプリカとのハミング距離を計算する。

［ステップＳ１２１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカのビットごとに、該当ビットをフリップした場合における遷移前後でのレプリカ間の相互作用のエネルギーの差分（ΔＧ_１，ΔＧ_２，・・・，ΔＧ_Ｎ）を計算する。例えば１番目のビットをフリップした場合のレプリカ間の相互作用のエネルギーの差分がΔＧ_１である。相互作用を与えるレプリカの１つ目の選択方法では、相互作用のエネルギーの差分は、ステップＳ１２０の処理で計算されたハミング距離を式（２７）に代入して得られる相互作用の強さを用いて計算される。

［ステップＳ１２２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカの規格化定数Ｚを計算する。例えばレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、レプリカ間距離がハミング距離の一次式の場合は、状態遷移前後での規格化定数の差分を計算してもよい。差分を計算した場合、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、状態遷移ごとの規格化定数の差分を積算することで、最新の規格化定数を得ることができる。

相互作用を与えるレプリカの１つ目の選択方法では、状態遷移を繰り返す処理が行われるたびに計算されるレプリカ間の距離の計算回数は、２ｓＭ回となる。そのため、ｓが小さければ、全レプリカ間の相互作用を考慮した場合の計算回数（Ｍ^２回）より大幅に計算回数を削減できる。

なお、全レプリカの状態を直接観測することはできないが、各レプリカに対して相互作用を与える範囲を±ｓのレプリカ番号の範囲とすることで、各レプリカ間の相互作用の影響が、レプリカ集団全体に波及することが期待できる。

（相互作用を与えるレプリカの２つ目の選択方法）
２つ目の選択方法は、各レプリカに与えられたラベルに基づいてＭ個のレプリカを複数のグループにグループ分けし、異なるグループ間に属するレプリカ間についてだけ、相互作用を与えるものである。この方法では、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカは、レプリカ番号＝ｌのレプリカが属するグループとは異なる各グループの代表レプリカに限られる。

２つ目の選択方法では、レプリカ番号＝ｌのレプリカに対して与えられる相互作用の強さは、式（１５）または式（１６）の代わりに、例えば、以下の式（２８）で定義できる。

式（２８）において、ｒはレプリカ番号＝ｌのレプリカが属するグループのグループ番号を表し、Ｒは全グループ数を表す。また、ｘ^（ｋ） _ｒｅｐは、グループ番号のｋのグループにおける代表レプリカを表す。式（２８）のように、レプリカ番号＝ｌのレプリカと、グループ番号＝ｒ以外のグループ番号のグループにおける代表レプリカとの距離に基づいて、相互作用の強さが定義される。

図１６は、相互作用を与えるレプリカの２つ目の選択方法の一例を示す図である。
図１６では、９個のレプリカが、３つのグループにグループ分けされている例が示されている。レプリカ番号＝１，２，３のレプリカは、グループ番号＝１のグループに属し、レプリカ番号＝４，５，６のレプリカは、グループ番号＝２のグループに属し、レプリカ番号＝７，８，９のレプリカは、グループ番号＝３のグループに属している。

また、図１６の例では、各グループにおいて、中間のレプリカ番号のレプリカが代表レプリカに設定されている。すなわち、グループ番号＝１のグループの代表レプリカは、レプリカ番号＝２のレプリカであり、グループ番号＝２のグループの代表レプリカは、レプリカ番号＝５のレプリカであり、グループ番号＝３のグループの代表レプリカは、レプリカ番号＝８のレプリカである。

図１６の例では、レプリカ番号＝１のレプリカに対して相互作用を与えるレプリカは、グループ番号＝２のグループに属するレプリカ番号＝５のレプリカと、グループ番号＝３のグループに属するレプリカ番号＝８のレプリカとなる。

図１７は、相互作用を与えるレプリカの２つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。以下、図１７に示す処理をステップ番号に沿って説明する。

［ステップＳ１３０］解探索エンジン３５０は、各レプリカがどのグループに属するかを示すグループ分け情報を、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。レプリカのグループ分けは、例えば、サーバ１００から全グループ数Ｒを与えられた制御装置２００によって予め行われており、各レプリカ番号に対して、グループ番号が紐付けられている。グループ分けの結果得られたグループ分け情報は、データ受け取り部３４０を介して解探索エンジン３５０に供給される。

［ステップＳ１３１］解探索エンジン３５０は、各グループの代表レプリカを示す情報をレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。各グループの代表レプリカは、例えば、制御装置２００によって決定される。例えば、図１６のように、各グループにおいて、真ん中のレプリカ番号のレプリカが代表レプリカとして決定される。決定された各グループの代表レプリカの情報は、データ受け取り部３４０を介して解探索エンジン３５０に供給される。

その後の処理（ステップＳ１３２からステップＳ１３９）は、ステップＳ１３３の処理以外、図１３の処理（ステップＳ１０１からステップＳ１０８）と同じである。
ステップＳ１３３のレプリカごとの解探索処理は、図１４に示した処理手順と同じ処理手順で行われるが、ステップＳ１１０の処理内容のうち、図１５に示したステップＳ１２０，Ｓ１２１の処理が、相互作用を与えるレプリカの１つ目の選択方法を適用した場合と異なる。例えば、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ステップＳ１２０，Ｓ１２１の処理において以下のような処理を行う。

ステップＳ１２０の処理において、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、グループ分け情報に基づいて、自身が解探索を担当するレプリカが属するグループを認識する。そして、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、代表レプリカを示す情報に基づいて、自身が解探索を担当するレプリカと、他の各グループの代表レプリカとの間のハミング距離を計算する。

ステップＳ１２１の処理において、相互作用を与えるレプリカの２つ目の選択方法では、相互作用のエネルギーの差分は、上記のように計算されたハミング距離を式（２８）に代入して得られる相互作用の強さを用いて計算される。

相互作用を与えるレプリカの２つ目の選択方法では、状態遷移を繰り返す処理が行われるたびに計算されるレプリカ間の距離の計算回数は、Ｍ（Ｒ－１）回となる。そのため、Ｒが小さければ、全レプリカ間の相互作用を考慮した場合の計算回数（Ｍ^２回）より大幅に計算回数を削減できる。

このような方法では、同じグループ内に属するレプリカの状態は同じように遷移し、グループごとに状態空間の探索が行われるようになる。
（相互作用を与えるレプリカの３つ目の選択方法）
３つ目の選択方法は、相互作用を適用するレプリカの範囲を動的に決定する方法である。この方法では、１つ目の方法と同じように、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカが、ｌ±ｓ_ｔのレプリカ番号の範囲のものとなるが、このｓ_ｔは動的に変化する。状態遷移を繰り返す処理が行われるたびに、レプリカ番号＝ｌのレプリカとｌ±ｓ_ｔのレプリカ番号の範囲の各レプリカとの距離の平均値が計算され、その平均値と２つの閾値（Ｄ_１，Ｄ_２（Ｄ_１＜Ｄ_２））との比較結果に基づいて、ｓ_ｔが減少または増加する。例えば、斥力の相互作用を発生させる場合、上記距離の平均値が、Ｄ_１よりも小さければｓ_ｔは＋１され、Ｄ_２よりも大きければｓ_ｔは－１される。引力の相互作用を発生させる場合はこの逆となる。なお、ｓ_ｔの変化幅は±１に限定されるわけではなく、±２でもよいし、それより大きい変化幅であってもよい。

３つ目の選択方法では、レプリカ番号＝ｌのレプリカに対して与えられる相互作用の強さは、式（１５）または式（１６）の代わりに、例えば、以下の式（２９）で定義できる。

式（２９）のように、レプリカ番号＝ｌのレプリカと、ｌ±ｓ_ｔのレプリカ番号の範囲のレプリカとの距離に基づいて、相互作用の強さが定義される。
なお、レプリカ番号＝ｌのレプリカとｌ±ｓ_ｔのレプリカ番号の範囲の各レプリカとの距離の平均値ｄ_ｔは、以下の式（３０）で表される。

図１８は、相互作用を与えるレプリカの３つ目の選択方法の一例を示す図である。
図１８では、１２個のレプリカに対して、１～１２のレプリカ番号が与えられており、このうち、各レプリカに相互作用を与えるレプリカの範囲（相互作用の適用範囲）を決めるｓ_ｔが、ある反復回数ｔであるときに、ｓ_ｔ＝２である例が示されている。

このとき、斥力の相互作用が強すぎる場合（上記の平均値ｄ_ｔ＞Ｄ_２の場合）、ｓ_ｔは－１され、次の反復回数ｔ＋１におけるｓ_ｔ＋１は、１となっている。
なお、１つ目の選択方法と同様に、負のレプリカ番号が発生することを避けるために、図１８の例のように、ｌ＝１の１つ前はｌ＝１２（ｌ＝１２の次はｌ＝１）というように、レプリカ番号が循環していると考える。つまり、ｋ＝ｌ±ｓ_ｔ ｍｏｄ Ｍが、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカの範囲であるものとする。

図１９は、相互作用を与えるレプリカの３つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。以下、図１９に示す処理をステップ番号に沿って説明する。

［ステップＳ１４０］解探索エンジン３５０は、前述の２つの閾値（Ｄ_１，Ｄ_２（Ｄ_１＜Ｄ_２））を、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。Ｄ_１，Ｄ_２は、例えば、サーバ１００によって決定され、制御装置２００に入力され、データ受け取り部３４０を介して解探索エンジン３５０に供給される。

［ステップＳ１４１］ステップＳ１４１の処理では、図１３のステップＳ１０１の処理と同様の処理が行われる。ただし、解探索エンジン３５０はさらに、ｔ＝０に初期化するとともに、上記のｓ_ｔの初期値ｓ_０をレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。初期値ｓ_０は、斥力の相互作用を発生させる場合はｓ_０＝０であり、引力の相互作用を発生させる場合はｓ_０＝１である。また、解探索エンジン３５０は上記の距離の平均値ｄ_ｔをｄ_０＝０に初期化する。

その後の処理（ステップＳ１４２からステップＳ１４８）のうち、ステップＳ１４２の処理以外、図１３の処理（ステップＳ１０３からステップＳ１０８）と同じである。
相互作用を与えるレプリカの３つ目の選択方法では、ステップＳ１４２のレプリカごとの解探索処理は、例えば、以下のように行われる。

図２０は、相互作用を与えるレプリカの３つ目の選択方法におけるレプリカごとの解探索処理の手順の一例を示すフローチャートである。以下、図２０に示す処理をステップ番号に沿って説明する。

［ステップＳ１５０］解探索エンジン３５０内のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｄ_ｔ＜Ｄ_１であるか否かを判断する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｄ_ｔ＜Ｄ_１であると判断した場合、処理をステップＳ１５１に進める。またレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｄ_ｔ＜Ｄ_１ではないと判断した場合、処理をステップＳ１５２に進める。

［ステップＳ１５１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｓ_ｔを更新する。相互作用として斥力を発生させる場合、ｓ_ｔ＝ｓ_ｔ＋１に更新され、相互作用として引力を発生させる場合、ｓ_ｔ＝ｓ_ｔ－１に更新される。

［ステップＳ１５２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｄ_ｔ＞Ｄ_２であるか否かを判断する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｄ_ｔ＞Ｄ_２であると判断した場合、処理をステップＳ１５３に進める。またレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｄ_ｔ＞Ｄ_２ではないと判断した場合、処理をステップＳ１５４に進める。

［ステップＳ１５３］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｓ_ｔを更新する。相互作用として斥力を発生させる場合、ｓ_ｔ＝ｓ_ｔ－１に更新され、相互作用として引力を発生させる場合、ｓ_ｔ＝ｓ_ｔ＋１に更新される。

その後の処理（ステップＳ１５４からステップＳ１６０）のうち、ステップＳ１５４の処理以外、図１４の処理（ステップＳ１１１からステップＳ１１６）と同じである。ステップＳ１６０の処理が終わると、ステップＳ１５６の処理でインクリメントされた反復回数ｔを用いて、ステップＳ１５０からの処理が繰り返される。

次に、ステップＳ１５４の処理である、レプリカ間の相互作用のエネルギーの差分（ΔＧ_１，ΔＧ_２，・・・，ΔＧ_Ｎ）の計算処理について詳細に説明する。
図２１は、レプリカ間の相互作用のエネルギーの差分の計算手順の一例を示すフローチャートである。以下、図２１に示す処理をステップ番号に沿って説明する。

［ステップＳ１７０］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、自身が解探索を担当するレプリカと、他の何れかのレプリカとの間のハミング距離を計算する。相互作用を与えるレプリカの３つ目の選択方法では、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、解探索を担当するレプリカと、そのレプリカのレプリカ番号に対して±ｓ_ｔの範囲のレプリカ番号のレプリカとのハミング距離を計算する。

［ステップＳ１７１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカのビットごとに、該当ビットをフリップした場合における遷移前後でのレプリカ間の相互作用のエネルギーの差分（ΔＧ_１，ΔＧ_２，・・・，ΔＧ_Ｎ）を計算する。相互作用を与えるレプリカの３つ目の選択方法では、相互作用のエネルギーの差分は、ステップＳ１７０の処理で計算されたハミング距離を式（２９）に代入して得られる相互作用の強さを用いて計算される。

［ステップＳ１７２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカの規格化定数Ｚを計算する。
［ステップＳ１７３］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ステップＳ１７０の処理で計算したハミング距離を用いて、式（３０）で表されるｄ_ｔを計算する。

相互作用を与えるレプリカの３つ目の選択方法では、状態遷移を繰り返す処理が行われるたびに計算されるレプリカ間の距離の計算回数は、２ｓ_ｔＭ回となる。そのため、ｓ_ｔが小さければ、全レプリカ間の相互作用を考慮した場合の計算回数（Ｍ^２回）より大幅に計算回数を削減できる。

斥力の相互作用を発生させる場合、レプリカ間の距離の平均が小さすぎる（ｄ_ｔ＜Ｄ_１）ということは、それらのレプリカは似たような状態にあり、意図した相互作用の効果を反映していないことの指標になる。そのため、次の反復回数においてはｓ_ｔが増加している。また、レプリカ間の距離の平均が大きすぎる（ｄ_ｔ＞Ｄ_２）ということは、それらのレプリカは大きく異なる状態にあり、これも意図した相互作用の効果を反映していないことの指標になる。そのため、次の反復回数においてはｓ_ｔが減少している。引力の相互作用を発生させる場合はこの反対である。このようにすることで、不必要な相互作用の発生を抑制できる。

（相互作用を与えるレプリカの４つ目の選択方法）
４つ目の方法は、相互作用を適用するレプリカの範囲をランダムに決定する方法である。この方法では、状態遷移を繰り返す処理が行われるたびに、レプリカ番号＝ｌのレプリカに対し、他の各レプリカについて、所定の確率ｐで相互作用を与えるレプリカとして採用する。ある反復回数ｔのときに、レプリカ番号＝ｌのレプリカに対して相互作用を与えるレプリカの範囲（レプリカの集合）をＣ_ｌ（ｔ）とすると、相互作用の強さは、式（１５）または式（１６）の代わりに、例えば、以下の式（３１）で定義できる。

図２２は、相互作用を与えるレプリカの４つ目の選択方法の一例を示す図である。
図２２では、４個のレプリカ間の相互作用の有無が示されている。図２２の例では、レプリカ番号＝１のレプリカに対して相互作用を与えるレプリカの範囲であるＣ_１（ｔ）は｛３，４｝である。すなわち、レプリカ番号＝３，４の２つのレプリカがレプリカ番号＝１のレプリカに対して相互作用を与える。また、レプリカ番号＝２のレプリカに対して相互作用を与えるレプリカの範囲であるＣ_２（ｔ）はφ（空集合を表す）である。すなわち、レプリカ番号＝２のレプリカに対して相互作用を与えるレプリカはない。また、レプリカ番号＝３のレプリカに対して相互作用を与えるレプリカの範囲であるＣ_３（ｔ）は｛１，４｝である。すなわち、レプリカ番号＝１，４の２つのレプリカがレプリカ番号＝３のレプリカに対して相互作用を与える。レプリカ番号＝４のレプリカに対して相互作用を与えるレプリカの範囲であるＣ_４（ｔ）は｛１，３｝である。すなわち、レプリカ番号＝１，３の２つのレプリカがレプリカ番号＝４のレプリカに対して相互作用を与える。

図２３は、相互作用を与えるレプリカの４つ目の選択方法を用いた場合の解探索処理の手順の一例を示すフローチャートである。以下、図２３に示す処理をステップ番号に沿って説明する。

［ステップＳ１８０］解探索エンジン３５０は、前述の確率ｐを、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。確率ｐは、例えば、サーバ１００によって決定され、制御装置２００に入力され、データ受け取り部３４０を介して解探索エンジン３５０に供給される。

［ステップＳ１８１］ステップＳ１８１の処理では、図１３のステップＳ１０１の処理と同様の処理が行われる。ただし、解探索エンジン３５０はさらに、ｔ＝０に初期化する。

その後の処理（ステップＳ１８２からステップＳ１８８）のうち、ステップＳ１８２の処理以外、図１３の処理（ステップＳ１０３からステップＳ１０８）と同じである。
相互作用を与えるレプリカの４つ目の選択方法では、ステップＳ１８２のレプリカごとの解探索処理は、例えば、以下のように行われる。

図２４は、相互作用を与えるレプリカの４つ目の選択方法におけるレプリカごとの解探索処理の手順の一例を示すフローチャートである。以下、図２４に示す処理をステップ番号に沿って説明する。

［ステップＳ１９０］解探索エンジン３５０内のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、前述のＣ_ｌ（ｔ）を計算する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、例えば、自身のレプリカのレプリカ番号をｌとした場合、Ｃ_ｌ（ｔ）＝φとする。そして、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ｌを除く各レプリカ番号に対して［０，１］の乱数Ｕ_１を与え、ｐ＜Ｕ_１ならばそのレプリカ番号をＣ_ｌ（ｔ）に加える処理をレプリカ番号の小さい順にレプリカ番号＝Ｍまで繰り返す。

その後の処理（ステップＳ１９１からステップＳ１９７）のうち、ステップＳ１９１の処理以外、図１４の処理（ステップＳ１１１からステップＳ１１６）と同じである。ステップＳ１９７の処理が終わると、ステップＳ１９３の処理でインクリメントされた反復回数ｔを用いて、ステップＳ１９０からの処理が繰り返される。

ステップＳ１９１のレプリカ間の相互作用のエネルギーの差分（ΔＧ₁，ΔＧ₂，・・・，ΔＧ_N）の計算処理の処理手順は図１５に示した処理手順と同じである。ただし、図１５のステップＳ１２０，Ｓ１２１の処理が、相互作用を与えるレプリカの１つ目の選択方法を適用した場合と異なる。例えば、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、ステップＳ１２０，Ｓ１２１の処理において以下のような処理を行う。

ステップＳ１２０の処理において、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、Ｃ_ｌ（ｔ）に含まれるレプリカ番号のレプリカと、自身が解探索を担当するレプリカ番号＝ｌのレプリカとの間のハミング距離を計算する。

ステップＳ１２１の処理において、相互作用を与えるレプリカの２つ目の選択方法では、相互作用のエネルギーの差分は、上記のように計算されたハミング距離を式（３１）に代入して得られる相互作用の強さを用いて計算される。

レプリカ間の距離の計算回数の期待値（平均の計算量）は、確率ｐを用いて以下の式（３２）のように表せる。

式（３２）において、ｉとｊはそれぞれレプリカ番号を表しており、“ｉ←→ｊ”は、レプリカ番号＝ｉ，ｊのレプリカ間に相互作用が与えられることを表している。Ｅは期待値、Ｐはレプリカ番号＝ｉ，ｊのレプリカ間に相互作用が与えられる確率を表している。また１｛ｉ←→ｊ｝はレプリカ番号＝ｉ，ｊのレプリカ間に相互作用が与えられる場合に１、相互作用が与えられない場合に０となる指示関数である。

式（３２）のように平均の計算量は、ｐＭ（Ｍ－１）／２となり、ｐの次数が１／Ｍとなる程度に小さければ、全レプリカ間の相互作用を考慮した場合の計算回数（Ｍ^２回）より大幅に計算回数を削減できる。

このような方法では、各反復回数において、相互作用が与えられる範囲をランダムに制限していることで、レプリカ番号の差が大きいレプリカ間であっても相互作用を与えられる可能性が高まる。つまり、レプリカ番号の差によらず、相互作用をレプリカ間に与える可能性があり、レプリカ番号による相互作用の適用範囲の偏りを抑制できる。

（更新ビット選択方法）
次に、図１４のステップＳ１１４、図２０のステップＳ１５８の更新ビットの選択方法について説明する。更新ビットの選択方法としては、例えば以下の３つの方法が考えられる。

第１の更新ビット選択方法は、Original Boltzmannの方法である。第２の更新ビット選択方法は、エネルギーの並列計算を行い、エネルギーが下がる方向を先に参照することで効率的にビット更新を行う方法である。第３の更新ビット選択方法は、１イテレーションで常にビットフリップが起きるようにしたRejection-freeの方法である。

図２５は、第１の更新ビット選択方法による更新ビット選択処理の手順の一例を示すフローチャートである。以下、図２５に示す処理をステップ番号に沿って説明する。
［ステップＳ２０１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、レプリカ間距離を考慮に入れた提案確率ｇ（ｘ^ｌ→ｘ^ｌ［ｊ］）に従ってビットｊを選択する。

［ステップＳ２０２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、メトロポリス基準の受け入れ確率ａ（ｘ^ｌ→ｘ^ｌ［ｊ］）に従って、選ばれたビットをフリップするか否かを判定する。

第１の更新ビット選択方法は単純な方法であり、計算が容易であるが、選択したビットのフリップの提案が棄却されることもある。提案が棄却された場合、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、図１４のステップＳ１１５（または図２０のステップＳ１５９）で「ＮＯ」と判断し、更新ビット選択処理を繰り返す。

第１の更新ビット選択方法は、提案分布に偏りがある影響を受けて受け入れ確率が小さくなってしまい、棄却ばかりが起こってしまう可能性がある。そこで更新ビットの提案が棄却された場合、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、オフセット値Ｅ_ｏｆｆの値を増加させることで、次回の更新ビットにおいて更新ビットが選択される確率を高めることができる。例えばレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、エネルギーが下がる方向が無くなる（エネルギー差がどのビット更新に対しても正になる）ときにはオフセット値Ｅ_ｏｆｆに所定の値を加算する。

またレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、エネルギーの並列計算を行い、エネルギーが下がる方向を先に参照することで効率的にビット更新を行う第２の更新ビット選択方法を適用することもできる。

図２６は、第２の更新ビット選択方法の処理手順の一例を示すフローチャートである。以下、図２６に示す処理をステップ番号に沿って説明する。
［ステップＳ２１１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、すべてのビットに対して、メトロポリス基準の受け入れ確率ａ（ｘ^ｌ→ｘ^ｌ［ｊ］）に従って、該当ビットが選択された場合にフリップするか否かを判定する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、各ビットに対応付けて、判定結果を示すフラグを設定する。

［ステップＳ２１２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、各ビットのフラグを参照し、ツリー状に接続されたセレクタを用いて、レプリカ間距離を考慮に入れた傾斜を与えて、更新ビットを選択する。

図２７は、更新ビット選択のためのツリー状に接続されたセレクタの一例を示す図である。制御回路３２０は、レプリカごとに、複数のビットそれぞれの状態遷移のエネルギーの変化値｛ΔＥ_ｉ｝に応じて、上記の式（６）、式（７）の受け入れ確率でその状態遷移を許容するか否かを判断する。そして、制御回路３２０は、状態遷移を受け入れると判断したビットのうちの１つを、ツリー状に接続されたセレクタによって選択する。制御回路３２０は、選択したビットの番号と、遷移可否Ｆとを出力する。

このように、制御回路３２０は、複数のビットそれぞれに対して並列探索を行うことで、更新ビットが選択できる確率を高めることができる。
並列探索を行うため、制御回路３２０は、次の回路構成を有する。一例として、ビットの数を３２個として説明する。図２７の例ではいずれか１つのビットのみが更新ビットとして選択されるものとする。

制御回路３２０は、比較回路部５１～５４とセレクタ部６０とを有する。
比較回路部５１～５４は、複数の状態変数のそれぞれが遷移した場合のエネルギーの変化値｛ΔＥ_ｉ｝を、ニューロン回路３１１，３１２，・・・，３１ｎから受け付ける。比較回路部５１～５４は、｛ΔＥ_ｉ｝に基づいて各状態遷移を受け入れるか否かを判定し、遷移可否｛ｆ_ｉ｝を出力する。比較回路部５１～５４それぞれは、８（＝３２／４）個の比較器を有する。比較回路部５１～５４に含まれる全ての比較器の合計数は３２個となる。

具体的には、比較回路部５１は、比較器Ｃ０，Ｃ１，…，Ｃ７を有する。比較回路部５２は、比較器Ｃ８，Ｃ９，…，Ｃ１５を有する。比較回路部５３は、比較器Ｃ１６，Ｃ１７，…，Ｃ２３を有する。比較回路部５４は、比較器Ｃ２４，Ｃ２５，…，Ｃ３１を有する。比較器Ｃｉ（図２７の例ではｉは０以上３１以下の整数）は、ΔＥ_ｉを受け付け、ΔＥ_ｉに基づく判定に応じて受け入れ可否ｆ_ｉを出力する。比較器Ｃｉによる判定では、エネルギーの変化値ΔＥ_ｉと温度パラメータＴの値を用いて算出した受け入れ確率と、乱数値ｕ（０≦ｕ≦１）とが比較される。例えば比較器Ｃｉは、乱数値ｕが受け入れ確率以下であれば、ｉ番目のビットのフリップを受け入れると判定する。

比較回路部５１～５４では、予め「Ｔ×ｌｏｇ（ｕ）」で表される値を計算することもできる。この値は、エネルギーが上がる状態遷移を確率的に引き起こす値であり、熱励起エネルギーまたは熱雑音と呼ぶこともできる。比較器Ｃｉは、ΔＥ_ｉと熱励起エネルギーとを比較し、例えば熱励起エネルギーの方が大きければ、ｉ番目のビットのフリップを受け入れると判定する。

セレクタ部６０には、比較器Ｃｉの出力値が状態遷移の候補として入力される。そしてセレクタ部６０は、複数の状態遷移の候補の何れか１つを選択し、出力する。セレクタ部６０は、当該選択を行うためのｎ段（ｎは２以上の整数）のセレクタツリーを有する。図２７の例では、ｎ＝５である。

セレクタツリーの第１段目（１ｓｔ）は、部分セレクタ部６０ａ，６０ｂを有する。セレクタツリーの第２段目（２ｎｄ）は、部分セレクタ部６０ｃを有する。セレクタツリーの第３段目（３ｒｄ）は、部分セレクタ部６０ｄを有する。セレクタツリーの第４段目（４ｔｈ）は、部分セレクタ部６０ｅを有する。セレクタツリーの第５段目（５ｔｈ）は、部分セレクタ部６０ｆを有する。

部分セレクタ部６０ａ，６０ｂ，…，６０ｆのそれぞれは、例えば、２つの入力のうちの１つを選択用乱数により選択して出力する１または複数のセレクタを有する。セレクタ６１は、複数のセレクタのうちの１つであり、他のセレクタもセレクタ６１と同様の構成である。セレクタ６１に対する２つの入力はｉとｊの遷移番号を特定するための識別値Ｎ_ｉ，Ｎ_ｊと遷移可否情報ｆ_ｉ，ｆ_ｊと提案確率ｇ（ｘ^ｌ→ｘ^ｌ［ｉ］），ｇ（ｘ^ｌ→ｘ^ｌ［ｊ］）である。セレクタ６１の出力は遷移可否情報ｆ_ｉ，ｆ_ｊの論理和として得られる可否情報ｆ_ｏと、ｉとｊのうち選択された方の遷移番号を特定するための識別値Ｎｏと、選択された方のビットの提案確率ｇ（ｘ^l→ｘ^l［ｏ］）である。

セレクタ６１は、遷移可否情報ｆ_ｉ，ｆ_ｊのいずれか一方が１（受け入れ可）、他方が０（受け入れ不可）の場合は受け入れ可の方のビットを必ず選択する。セレクタ６１は、遷移可否情報ｆ_ｉ，ｆ_ｊの両方０の場合はどのように選んでもよい。

セレクタ６１は、遷移可否情報ｆ_ｉ，ｆ_ｊの両方が１の場合には、候補選択用乱数を用いて、提案確率に応じた確率で一方を選択する。例えばセレクタ６１は、提案確率ｇ（ｘ^ｌ→ｘ^ｌ［ｉ］），ｇ（ｘ^ｌ→ｘ^ｌ［ｊ］）の比率に応じて、０から１の値域を、ｉとｊのビットに対応する２つの区間に分ける。そしてセレクタ６１は、候補選択用乱数を含む区間に対応するビットを選択する。そしてセレクタ６１は、選択結果により選ばれたビットの識別値Ｎ_ｏを生成し出力する。

図２７の例ではセレクタ６１以外のセレクタが略記されている。図２７では黒い丸印で表された箇所が、１つのセレクタに相当する。部分セレクタ部６０ａ，６０ｂ，６０ｃのそれぞれは、８個のセレクタを有する。部分セレクタ部６０ｄは、４個のセレクタを有する。部分セレクタ部６０ｅは、２個のセレクタを有する。部分セレクタ部６０ｆは、１個のセレクタを有する。部分セレクタ部６０ａ～６０ｆ内の各セレクタがセレクタ６１と同様の選択処理を行うことで、レプリカ間距離に応じた提案確率が高いビットほど選択される可能性を高くして、１つのビットが状態遷移の候補として出力される。

図２７で示されるように、制御回路３２０は、状態遷移の並列探索を行い、セレクタの２進木構造を用いてノックダウン方式（あるいはトーナメント方式とも呼ばれる）で、状態遷移の候補を１つに絞り込む。フリップによりエネルギーが減少するビットは、比較器により受け入れ可と判断されるため、フリップによりエネルギーが減少するビットが少なくとも１つ存在すれば、セレクタ部６０による１回の選択で更新ビットを選択できる。また局所解に達しており、いずれのビットをフリップしてもエネルギーが増加する場合であっても、乱数値ｕと温度パラメータＴの値に基づいて、いずれか１つのビットのフリップが受け入れられる可能性がある。いずれか１つのビットのフリップが受け入れられれば、セレクタ部６０による１回の選択で更新ビットを選択できる。しかも、セレクタによる選択時に、レプリカ間距離を反映させた提案確率を用いたことで、提案確率が高いビットほど更新ビットとして選択される可能性が高くなる。

なおレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、セレクタ部６０が出力した遷移可否情報が０（受け入れ不可）の場合、オフセット値を増加させて、更新ビット選択処理を繰り返す。これにより、更新ビットを早期に選択できる可能性が高くなる。

図２８は、第３の更新ビット選択方法の処理手順の一例を示すフローチャートである。第３の更新ビット選択方法は、以下の１ステップで更新ビットを選択できる。
［ステップＳ２３１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、各ビットの遷移確率Ｗ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）＝ｇ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）×ａ（ｘ^ｌ→ｘ^ｌ［ｊ_０］を用いて、以下の式（３３）に示すRejection-freeの遷移確率Ｗ（ｘ^ｌ→ｘ^ｌ［ｊ_０］）（Ｗはチルダ付き）を計算する。

そしてレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、Ｒｅｊｅｃｔｉｏｎ－ｆｒｅｅの遷移確率により、いずれか１つのビットを更新ビットとして選択する。このように各ビットの遷移確率を正規化し、受け入れ確率の合計が１となるようにすることで、１回の更新ビット選択処理で、常に更新ビットを選択することが可能となる。

以上説明したように、第２の実施の形態に係るイジングマシン３００は、レプリカ間の相互作用を提案確率に反映させ、複数のレプリカを用いた解探索を行っている。これにより、組み合わせ最適化問題をメトロポリスヘイスティングの方法に基づいて求解する際に、収束先の分布を保ったまま、それぞれのレプリカがバラバラに状態空間を探索することが期待され、求解性能が向上する。すなわち、最適解へ到達する可能性が高くなり、エネルギーの下がり方を速くすることができる。

また、イジングマシン３００は、全レプリカ間の相互作用を考慮するのではなく、一部のレプリカ間の相互作用を考慮するものであるため、全レプリカ間の相互作用を考慮した場合のレプリカ間の距離の計算回数（Ｍ^２回）よりも計算回数を削減できる。例えば、前述の相互作用を与えるレプリカの４つの選択方法によれば、上記計算回数をＭの１次式で表せる程度に抑えることができる。

図２９は、レプリカ間に斥力の相互作用を設定した場合のエネルギーランドスケープを示す図である。複数のレプリカ７１～７３のうち、レプリカ７１，７２の間と、レプリカ７２，７３の間には、斥力の相互作用が与えられている。レプリカ７１，７３の間には、斥力の相互作用が与えられていない。この場合、レプリカ７１，７２とレプリカ７２，７３が互いに反発し合うことで、広い探索空間を効率的に探索することができる。

図３０は、レプリカ間に引力の相互作用を設定した場合のエネルギーランドスケープを示す図である。複数のレプリカ７４，７５の間には、引力の相互作用が与えられている。レプリカ７４，７５が互いに引きつけられることで局所解から脱出し易くなり、集団全体として大域解に到達できる可能性が高まる。一方、レプリカ７４，７６間には引力の相互作用が与えられていない。この場合、レプリカ７６がレプリカ７４に引き寄せられて局所解に嵌まることが抑制される。

次に、効果確認を行った検証例について、図３１～図３３を参照して説明する。
図３１は、第１の検証例を示す図である。図３２および図３３は、第２の検証例を示す図である。

図３１～図３３に示す例は、二次割り当て問題（ＱＡＰ）という代表的な組み合わせ最適化問題のいくつかのインスタンスについて検証した結果である。提案分布に従った各ビットの提案確率の計算には前述の式（１７）を使用している。レプリカ間の相互作用のエネルギーとしては、前述の式（１９）に示したハミング距離の一次関数を使用している。更新ビットの選択方法としては、第３の更新ビット選択方法（Rejection-free）が用いられている。また、全レプリカ数Ｍ＝３０である。

図３１の例では、１－ｂｉｔフリップ遷移かつレプリカ交換を使う解探索手法において、レプリカ間の相互作用の有無によるエネルギーの下がり方の違いを比較している。横軸が状態遷移の反復回数であり、縦軸がその時点で得られているエネルギーの最小値である。γ（図３１では“ｇａｍｍａ”と表記）を斥力相互作用のパラメータとしたときに、ｇａｍｍａ＝０とｇａｍｍａ＜０（つまり斥力相互作用の有無）との場合についてエネルギーの下がり方について比較している。

図３１の例では、斥力相互作用を導入した場合（ｇａｍｍａ－３）のほうが相互作用を導入しない場合（ｇａｍｍａ－０）よりもエネルギーの下がり方が速い。
このように、レプリカ間の相互作用を導入したことで、解探索性能が向上している。しかもレプリカ間の相互作用を提案確率に反映しており、目的関数に手を加えないため、適切な目的関数（例えばギブス分布を示す関数）を用いた解探索が可能となる。

図３２および図３３に示す例は、相互作用を与えるレプリカの１つ目の選択方法を用いたものであり、相互作用の適用範囲を決める前述のｓを、１～１５の範囲で変えた場合について、エネルギーの下がり方の違いを比較したものである。なお、ｓ＝０は、相互作用を適用しない場合を示している。

図３２、図３３に示すように、相互作用の適用範囲が広いほう（図３３）が、エネルギーの下がり方が速いが、相互作用の適用範囲が狭い場合（図３２）においても、相互作用を適用しない場合よりも低いエネルギーへの収束が生じている。

なお、非特許文献１に示されたＣＭＣと呼ばれる方法は、実数を定義域とする目的関数にのみ適用可能な方法であり、２値の離散空間を定義域とする（バイナリ変数）イジングマシンの目的関数に直接は適応できない。またＣＭＣでは、距離が近いレプリカの数（密度）をカウントしているが、１－ｂｉｔフリップの場合にレプリカすべての状態をみたときに、その状況がフリップした前後で大きくは変わらない。そのため、あるビットのフリップ前後でのレプリカ数の密度の比はほぼ１に近くなってしまい、２値の離散空間を定義域とするとレプリカ相互作用の効果が薄くなってしまう。それに対して、第２の実施の形態に示した方法では、２値の離散空間を定義域とする組み合わせ最適化問題に適用でき、求解性能も向上する。

また、非特許文献２に示されたＲＥと呼ばれる方法では、レプリカ間の相互作用を目的関数に直接加える方法をとっているため、本来の目的関数の最適化を行っている保証はない。それに対して、第２の実施の形態に示した方法では、レプリカ間の相互作用を提案確率に反映しており、適切な目的関数を用いた解探索が可能となる。

〔その他の実施の形態〕
第２の実施の形態では、レプリカ間の温度交換を行っているが、レプリカ間の温度交換を実施せずに、複数のレプリカで個別に解探索を行うことも可能である。その場合であっても、レプリカ間の相互作用を考慮した解探索により、解探索能力が向上する。

また第２の実施の形態では、２値の離散空間を定義域とするイジングモデルを用いた求解を行っているが、実数を定義域とするモデルをレプリカとして求解する場合にも適用可能である。

さらに第２の実施の形態では、複数のニューロン回路３１１，３１２，・・・，３１ｎを有するイジングマシン３００で解探索を行っているが、同じ処理を図３に示したサーバ１００と同様のハードウェア構成のノイマン型コンピュータで実現することも可能である。その場合、イジングマシン３００は、例えばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態と同様の解探索処理を実行する。イジングマシン３００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。例えば、イジングマシン３００に実行させるプログラムをストレージ装置に格納しておくことができる。イジングマシン３００のプロセッサは、ストレージ装置内のプログラムの少なくとも一部をメモリにロードし、プログラムを実行する。またイジングマシン３００に実行させるプログラムを、光ディスク、メモリ装置、メモリカードなどの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、例えばイジングマシン３００のプロセッサからの制御により、ストレージ装置にインストールされた後、実行可能となる。またイジングマシン３００のプロセッサが、可搬型記録媒体から直接プログラムを読み出して実行することもできる。

以上、実施の形態を例示したが、実施の形態で示した各部の構成は同様の機能を有する他のものに置換することができる。また、他の任意の構成物や工程が付加されてもよい。さらに、前述した実施の形態のうちの任意の２以上の構成（特徴）を組み合わせたものであってもよい。

１ 状態空間
２～４ レプリカ
１０ 最適化装置
１１ 記憶部
１２ 処理部

Claims (10)

1. 複数のレプリカそれぞれの複数の状態変数の値を記憶する記憶部と、
   前記複数のレプリカそれぞれについて、該レプリカが有する前記複数の状態変数のうちの第１状態変数の値を更新した場合における、前記レプリカと、前記複数のレプリカから前記レプリカを除いたレプリカ群のうちの一部である他のレプリカとの、前記複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での距離の変化に応じた相互作用の強さの変化量を特定し、前記第１状態変数の値を更新した場合における前記相互作用の強さの前記変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、前記第１状態変数の値を更新するか否かを決定する処理部と、
   を有する最適化装置。
2. 前記複数のレプリカそれぞれには、前記複数のレプリカそれぞれを識別する識別情報であるレプリカ番号が与えられており、
   前記他のレプリカは、前記レプリカを識別する識別情報である第１のレプリカ番号との差が所定の値である範囲に含まれる第２のレプリカ番号が与えられているレプリカである、
   請求項１に記載の最適化装置。
3. 前記処理部は、前記距離の平均値と、第１の閾値または前記第１の閾値より大きい第２の閾値との比較結果に基づいて、前記所定の値を変化させる、請求項２に記載の最適化装置。
4. 前記処理部は、斥力の前記相互作用を発生させる場合、前記平均値が前記第１の閾値よりも小さいときは、前記所定の値を増加させ、前記平均値が前記第２の閾値よりも大きいときは、前記所定の値を減少させる、請求項３に記載の最適化装置。
5. 前記処理部は、引力の前記相互作用を発生させる場合、前記平均値が前記第１の閾値よりも小さいときは、前記所定の値を減少させ、前記平均値が前記第２の閾値よりも大きいときは、前記所定の値を増加させる、請求項３に記載の最適化装置。
6. 前記所定の値をｓ、前記複数のレプリカの数をＭ、前記第１のレプリカ番号をｌとし、前記レプリカ番号が循環しているとしたとき、前記他のレプリカは、ｌ－ｓをＭで割った余りと、ｌ＋ｓをＭで割った余りの範囲に含まれる前記第２のレプリカ番号が与えられているレプリカである請求項２ないし５の何れか一項に記載の最適化装置。
7. 前記複数のレプリカは、複数のグループにグループ分けされており、前記複数のグループそれぞれには代表レプリカが設定されており、
   前記複数のグループのうち第１のグループに属する前記レプリカに対する前記他のレプリカは、前記複数のグループのうち前記第１のグループ以外の他のグループの前記代表レプリカである、
   請求項１に記載の最適化装置。
8. 前記処理部は、前記レプリカ以外の前記複数のレプリカそれぞれについて、所定の確率で前記他のレプリカとして採用する、請求項１に記載の最適化装置。
9. 最適化装置が、
   複数の状態変数を有する複数のレプリカそれぞれについて、該レプリカが有する前記複数の状態変数のうちの第１状態変数の値を更新した場合における、前記レプリカと、前記複数のレプリカから前記レプリカを除いたレプリカ群のうちの一部である他のレプリカとの、前記複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での距離の変化に応じた相互作用の強さの変化量を特定し、前記第１状態変数の値を更新した場合における前記相互作用の強さの前記変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、前記第１状態変数の値を更新するか否かを決定する、
   最適化方法。
10. 最適化装置に、
    複数の状態変数を有する複数のレプリカそれぞれについて、該レプリカが有する前記複数の状態変数のうちの第１状態変数の値を更新した場合における、前記レプリカと、前記複数のレプリカから前記レプリカを除いたレプリカ群のうちの一部である他のレプリカとの、前記複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での距離の変化に応じた相互作用の強さの変化量を計算し、前記第１状態変数の値を更新した場合における前記相互作用の強さの前記変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、前記第１状態変数の値を更新するか否かを決定する、
    処理を実行させる最適化プログラム。

Images