(57)【要約】【課題】計算処理の高速性を損なうことなく高精度のサンプリングを行う。【解決手段】オフセット制御部１３が、エネルギー変化ΔＥ１〜ΔＥＮに基づいて各状態遷移のうちの少なくとも１つが許容されるようにオフセット値Ｅｏｆｆを決定し、エネルギー変化ΔＥ１〜ΔＥＮのそれぞれにオフセット値Ｅｏｆｆを付加した複数の評価値を出力するとともに、所定試行回数毎にオフセット値Ｅｏｆｆを出力し、比較部１４が、温度値Ｔと乱数値とに基づいて決まる閾値と、複数の評価値のそれぞれとを比較した結果に基づき、各状態遷移について状態遷移を許容するか否かを示すフラグ値Ｆ１〜ＦＮを出力し、フラグ計数部１５が、フラグ値Ｆ１〜ＦＮのうち、状態遷移を許容することを示すフラグ値の数を計数した計数値Ｎｆを、所定試行回数毎に出力する。【選択図】図１
2020135727.tif 000002
【特許請求の範囲】
【請求項１】
エネルギーを表す評価関数に含まれる複数の状態変数の値をそれぞれ保持するとともに、所定試行回数毎に前記複数の状態変数の値を出力する状態保持部と、
前記複数の状態変数の値の何れかが変化することに応じて状態遷移が起こる場合、更新インデックス値に基づいて選択された重み値に基づき、エネルギーの変化値を各状態遷移に対して計算するエネルギー変化計算部と、
複数の前記エネルギーの変化値に基づいて前記各状態遷移のうちの少なくとも１つが許容されるようにオフセット値を決定し、複数の前記エネルギーの変化値のそれぞれに前記オフセット値を付加した複数の評価値を出力するとともに、前記所定試行回数毎に前記オフセット値を出力するオフセット制御部と、
温度を示す温度値と乱数値とに基づいて決まる閾値と、前記複数の評価値のそれぞれとを比較した結果に基づき、前記各状態遷移について状態遷移を許容するか否かを示す複数のフラグ値を出力する比較部と、
前記複数のフラグ値のうち、状態遷移を許容することを示すフラグ値の数を計数した計数値を、前記所定試行回数毎に出力するフラグ計数部と、
前記複数のフラグ値、または複数の前記エネルギーの変化値に基づいて、前記各状態遷移のうちの１つに対応するインデックス値を前記更新インデックス値として出力する選択部と、
を有するサンプリング装置。
【請求項２】
前記オフセット制御部が出力した前記オフセット値と、前記フラグ計数部が出力した前記計数値に基づき、マルコフ連鎖モンテカルロ法にしたがった状態遷移における１つの状態に滞在する試行回数の期待値を計算する試行回数計算部を、さらに有する請求項１に記載のサンプリング装置。
【請求項３】
前記オフセット制御部は、複数の前記エネルギーの変化値のうちの最小値が０以上である場合、前記最小値を前記オフセット値として決定し、複数の前記エネルギーの変化値のそれぞれから前記最小値を差し引くことで、前記複数の評価値を算出する、請求項１または２に記載のサンプリング装置。
【請求項４】
前記選択部は、前記複数のフラグ値のそれぞれに、それぞれ独立な複数の乱数値の何れかを加算した加算結果の大小関係に基づき、前記複数のフラグ値の何れかを選択する、請求項１乃至３の何れか一項に記載のサンプリング装置。
【請求項５】
前記複数のフラグ値に基づき、複数の前記エネルギーの変化値のうち、許容されない状態遷移であるとされている複数の状態遷移に対する複数の非許容エネルギー変化値を出力するエネルギー変化置き換え部と、
前記複数の非許容エネルギー変化値に基づき前記複数の状態遷移のうちの少なくとも１つが許容されるように前記オフセット値とは別のオフセット値を決定し、前記複数の非許容エネルギー変化値のそれぞれに前記別のオフセット値を付加した、前記複数の評価値とは別の複数の評価値を出力するとともに、前記所定試行回数毎に前記別のオフセット値を出力する他のオフセット制御部と、
前記閾値と、前記別の複数の評価値のそれぞれとを比較した結果に基づき、前記複数の状態遷移のそれぞれについて状態遷移を許容するか否かを示す、前記複数のフラグ値とは別の複数のフラグ値を出力する他の比較部と、
前記別の複数のフラグ値のうち、状態遷移を許容することを示すフラグ値の数を計数した、前記計数値とは別の計数値を、前記所定試行回数毎に出力する他のフラグ計数部と、
前記別の計数値、前記オフセット値及び前記別のオフセット値に基づいて、前記計数値を補正する計数値補正部と、
をさらに有する請求項１に記載のサンプリング装置。
【請求項６】
補正された前記計数値と、前記オフセット値と、前記別のオフセット値と、に基づいて、マルコフ連鎖モンテカルロ法にしたがった状態遷移における１つの状態に滞在する試行回数の期待値を計算する試行回数計算部を、さらに有する請求項５に記載のサンプリング装置。
【請求項７】
前記選択部は、複数の前記エネルギーの変化値のうち、０より小さい変化値を０に更新させた複数の更新エネルギー変化値を生成し、前記複数の更新エネルギー変化値のそれぞれに、前記温度値と互いに独立な複数の乱数値とに基づいてそれぞれ算出される複数の第２の閾値の何れかを加算した複数の第２の評価値のうちで、最小となる第２の評価値に対応する前記インデックス値を前記更新インデックス値として出力する、請求項１乃至６の何れか一項に記載のサンプリング装置。
【請求項８】
サンプリング装置の制御方法において、
前記サンプリング装置が有する状態保持部が、エネルギーを表す評価関数に含まれる複数の状態変数の値をそれぞれ保持するとともに、所定試行回数毎に前記複数の状態変数の値を出力し、
前記サンプリング装置が有するエネルギー変化計算部が、前記複数の状態変数の値の何れかが変化することに応じて状態遷移が起こる場合、更新インデックス値に基づいて選択された重み値に基づき、エネルギーの変化値を各状態遷移に対して計算し、
前記サンプリング装置が有するオフセット制御部が、複数の前記エネルギーの変化値に基づいて前記各状態遷移のうちの少なくとも１つが許容されるようにオフセット値を決定し、複数の前記エネルギーの変化値のそれぞれに前記オフセット値を付加した複数の評価値を出力するとともに、前記所定試行回数毎に前記オフセット値を出力し、
前記サンプリング装置が有する比較部が、温度を示す温度値と乱数値とに基づいて決まる閾値と、前記複数の評価値のそれぞれとを比較した結果に基づき、前記各状態遷移について状態遷移を許容するか否かを示す複数のフラグ値を出力し、
前記サンプリング装置が有するフラグ計数部が、前記複数のフラグ値のうち、状態遷移を許容することを示すフラグ値の数を計数した計数値を、前記所定試行回数毎に出力し、
前記サンプリング装置が有する選択部が、前記複数のフラグ値、または複数の前記エネルギーの変化値に基づいて、前記各状態遷移のうちの１つに対応するインデックス値を前記更新インデックス値として出力する、
サンプリング装置の制御方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、サンプリング装置及びサンプリング装置の制御方法に関する。
【背景技術】
【０００２】
組合せ最適化問題を解くための手法として、組合せ最適化問題を磁性体のスピンの振る舞いを表すイジングモデルに変換し、マルコフ連鎖モンテカルロ法を用いて、イジングモデルの状態をエネルギーが低い状態に遷移させていく手法がある。以下、マルコフ連鎖モンテカルロ法を、ＭＣＭＣ（Markov-Chain Monte Carlo）法と略す。イジングモデルの状態は、複数の状態変数の値の組合せにより表現でき、状態変数の数がＮ個の場合、ｘ＝（ｘ１，ｘ２，…，ｘＮ）などと表せる。各状態変数の値として、０または１を用いることができる。
【０００３】
イジングモデルのエネルギーを表すイジング型のエネルギー関数Ｅ（ｘ）は、たとえば、以下の式（１）で定義される。
【０００４】
【数１】
2020135727.tif 000003
【０００５】
右辺の１項目は、イジングモデルの全状態変数の全組合せについて、漏れと重複なく、２つの状態変数の値（０または１）と重み値（２つの状態変数の間の相互作用の強さを表す）との積を積算したものである。ｘｉは、識別情報（以下インデックス値という）がｉの状態変数、ｘｊは、インデックス値＝ｊの状態変数であり、Ｗｉｊは、インデックス値＝ｉ，ｊの状態変数間の相互作用の大きさを示す重み値である。
【０００６】
右辺の２項目は、各インデックス値についてのバイアス係数と状態変数との積の総和を求めたものである。ｂｉは、インデックス値＝ｉについてのバイアス係数を示している。
また、状態変数ｘｉが変化して１−ｘｉとなると、状態変数ｘｉの増加分は、Δｘｉ＝（１−ｘｉ）−ｘｉ＝１−２ｘｉと表せる。スピン反転（状態変数の値の変化）に伴うエネルギー変化ΔＥｉは、以下の式（２）で表される。
【０００７】
【数２】
2020135727.tif 000004
【０００８】
式（２）において、状態変数ｘｉが１から０に変化するとき、Δｘｉは−１となり、状態変数ｘｉが０から１に変化するとき、Δｘｉは１となる。なお、ｈｉはローカルフィールド（局所場）と呼ばれ、Δｘｉに応じてローカルフィールドｈｉに符号（＋１または−１）を乗じたものがエネルギー変化ΔＥｉとなる。
【０００９】
ある状態遷移に伴うエネルギー変化ΔＥに対するその状態遷移の許容確率として、以下の式（３）で表されるメトロポリス法またはギブス法で規定される許容確率Ａ（ΔＥ）を用いることができる。
【００１０】
【数３】
2020135727.tif 000005
【００１１】
式（３）においてβは、逆温度（温度を表す温度値の逆数）である。エネルギーが増加する状態遷移についても確率的に許容される。
通常のＭＣＭＣ法では、ランダムまたはインデックス値の順（シーケンシャル）に、状態変数が選択され、その状態変数の値が変化する状態遷移に伴うエネルギー変化ΔＥに基づいて、上記の許容確率Ａ（ΔＥ）でその状態遷移を許容する。そしてその状態遷移を許容する場合には状態変数の値が更新される。このような処理が所定の試行回数、繰り返される。また、最低エネルギーとなる状態（最適解）を探索するために、温度を徐々に下げていくシミュレーテッド・アニーリング法が用いられる場合もある。
【００１２】
一方、組合せ最適化問題を高速に解くためのハードウェアとして、デジタル回路を用いた最適化装置がある（たとえば、特許文献１参照）。最適化装置は、以下に示すように複数の状態遷移を同時に遷移候補として、１つの状態遷移を選択する並列探索を行う。
【００１３】
最適化装置は、各状態遷移に伴うエネルギー変化と温度値に基づいて、上記の許容確率Ａ（ΔＥ）で各状態遷移を許容する。そして、最適化装置は複数の状態遷移の中から、許容された状態遷移を優先的に１つ選択し状態を更新する。最適化装置は、以上のような処理を、所定の試行回数、繰り返し行う。また、従来の最適化装置では、状態が局所解に拘束されるなどによって状態遷移が生じない場合、エネルギー変化にオフセット値を付加していた。
【先行技術文献】
【特許文献】
【００１４】
【特許文献１】特開２０１８−６３６２６号公報
【発明の概要】
【発明が解決しようとする課題】
【００１５】
ところで、平衡状態では各状態の占有確率を示す確率分布はボルツマン分布となる。そのため、ＭＣＭＣ法を用いて温度一定とした状態遷移を繰り返す過程により得られた状態または状態に基づいた値をサンプルとして出力することで、ボルツマン分布にしたがうサンプルを発生するサンプラーが実現できる。発生したサンプルは、たとえば、機械学習などでの期待値の計算に用いられる。
【００１６】
しかし、従来のデジタル回路を用いた最適化装置は、前述のような並列探索処理及びオフセット値の付加による計算処理の高速化を行っているため、確率分布がボルツマン分布からずれる可能性がある。従来のデジタル回路を用いた最適化装置をサンプラー（以下サンプリング装置という）として用いる場合、確率分布のボルツマン分布からのずれによって、サンプリング精度が通常のＭＣＭＣ法を用いた場合よりも悪化してしまう可能性があった。
【００１７】
１つの側面では、本発明は、計算処理の高速性を損なうことなく高精度のサンプリングが可能なサンプリング装置及びサンプリング装置の制御方法を提供することを目的とする。
【課題を解決するための手段】
【００１８】
１つの実施態様では、エネルギーを表す評価関数に含まれる複数の状態変数の値をそれぞれ保持するとともに、所定試行回数毎に前記複数の状態変数の値を出力する状態保持部と、前記複数の状態変数の値の何れかが変化することに応じて状態遷移が起こる場合、更新インデックス値に基づいて選択された重み値に基づき、エネルギーの変化値を各状態遷移に対して計算するエネルギー変化計算部と、複数の前記エネルギーの変化値に基づいて前記各状態遷移のうちの少なくとも１つが許容されるようにオフセット値を決定し、複数の前記エネルギーの変化値のそれぞれに前記オフセット値を付加した複数の評価値を出力するとともに、前記所定試行回数毎に前記オフセット値を出力するオフセット制御部と、温度を示す温度値と乱数値とに基づいて決まる閾値と、前記複数の評価値のそれぞれとを比較した結果に基づき、前記各状態遷移について状態遷移を許容するか否かを示す複数のフラグ値を出力する比較部と、前記複数のフラグ値のうち、状態遷移を許容することを示すフラグ値の数を計数した計数値を、前記所定試行回数毎に出力するフラグ計数部と、前記複数のフラグ値、または複数の前記エネルギーの変化値に基づいて、前記各状態遷移のうちの１つに対応するインデックス値を前記更新インデックス値として出力する選択部と、を有するサンプリング装置が提供される。
【００１９】
また、１つの実施態様では、サンプリング装置の制御方法が提供される。
【発明の効果】
【００２０】
１つの側面では、本発明は、計算処理の高速性を損なうことなく高精度のサンプリングが可能となる。
【図面の簡単な説明】
【００２１】
【図１】第１の実施の形態のサンプリング装置の一例を示す図である。
【図２】各状態に留まる試行回数を１回とした場合の選択確率の分布例を示す図である。
【図３】状態保持部及びエネルギー変化計算部の一例を示す図である。
【図４】オフセット制御部及び比較部の一例を示す図である。
【図５】オフセット値の適用例を示す図である。
【図６】フラグ計数部の一例を示す図である。
【図７】フラグ計数部の他の例を示す図である。
【図８】選択部の一例を示す図である。
【図９】サンプリングタイミングの例を示す図である。
【図１０】サンプリング動作の一例の流れを示すフローチャートである。
【図１１】最適化処理動作の一例の流れを示すフローチャートである。
【図１２】第２の実施の形態のサンプリング装置の一例を示す図である。
【図１３】２種類のフラグ値についての計数値Ｎｆ１，Ｎｆ２を発生させる動作の一例の流れを示すフローチャートである。
【図１４】第３の実施の形態のサンプリング装置の一例を示す図である。
【図１５】閾値生成部の一例を示す図である。
【図１６】温度とエネルギーとの関係についてのシミュレーション結果の例を示す図である。
【図１７】試行回数と最低エネルギー状態が得られる確率との関係についてのシミュレーション結果の例を示す図である。
【図１８】レプリカ交換法を利用したサンプリング装置の例を示す図である。
【図１９】レプリカ交換法を利用したサンプリング装置の動作例を示すフローチャートである。
【図２０】サンプリング装置の他の変形例を示す図である。
【図２１】期待値計算処理を行うサンプリング装置の動作例を示すフローチャートである。
【発明を実施するための形態】
【００２２】
以下、発明を実施するための形態を、図面を参照しつつ説明する。
なお、イジングモデルを用いたシミュレーテッド・アニーリングにおいては、状態遷移に伴い変化する状態変数は１つだけである。そこで、以下では各状態遷移をそれぞれ識別するインデックス値は、１つの状態変数のインデックス値と等しいものとして説明を行う。しかし、状態遷移のインデックス値と状態遷移に伴い変化する状態遷移のインデックス値が一致する形態に限定されるものではない。
【００２３】
従来の最適化装置において確率分布がボルツマン分布からずれる原因は、並列探索処理及びオフセット値の付加による計算処理の高速化を行っているため、１つの状態に滞在する試行回数が、通常のＭＣＭＣ法を用いた場合よりも少なくなるためである。従来の最適化装置において、理想的なオフセット値が適用された場合、１つの状態に滞在する試行回数は、１回となる。
【００２４】
ボルツマン分布にしたがったサンプルを得るために、上記のようなオフセット値を用い、サンプリングを行うたび、通常のＭＣＭＣ法であればその状態に留まったであろう試行回数の値またはその期待値により、サンプルに重み付けを行えばよい。つまり、ある状態が１回観測された（サンプリングされた）場合、その状態が通常のＭＣＭＣ法を用いた場合の試行回数だけサンプリングされたとするような重み付けを行えばよい。
【００２５】
通常のＭＣＭＣ法において、ある状態からインデックス値＝ｉの状態変数の値が変化した別の状態に遷移する確率は、Γｉ＝Ａ（ΔＥｉ）／Ｎである。Ａ（ΔＥｉ）は、インデックス値＝ｉの状態変数の値が変化する状態遷移の許容確率である。Ｎは、状態変数の総数である。試行毎に別の状態に遷移する確率は、以下の式（４）で表せる。
【００２６】
【数４】
2020135727.tif 000006
【００２７】
式（４）で表される確率の逆数が、通常のＭＣＭＣ法において、１つの状態に滞在する試行回数の期待値となる。
なお、何回かの試行を行って最終的に状態変数ｘｉの値の変化が生じる確率Ｐｉは、以下の式（５）で表せる。
【００２８】
【数５】
2020135727.tif 000007
【００２９】
以下に示す第１の実施の形態のサンプリング装置は、所定試行回数毎に、通常のＭＣＭＣ法において１つの状態に滞在する試行回数の期待値を、その状態とともに出力する装置である。
【００３０】
（第１の実施の形態）
図１は、第１の実施の形態のサンプリング装置の一例を示す図である。
サンプリング装置１０は、状態保持部１１、エネルギー変化計算部１２、オフセット制御部１３、比較部１４、フラグ計数部１５、選択部１６、制御部１７、試行回数計算部１８を有する。
【００３１】
状態保持部１１は、たとえば、エネルギーを表す評価関数に含まれる複数の状態変数の値をそれぞれ保持するとともに、所定試行回数毎に複数の状態変数の値を出力する。評価関数は、たとえば、式（１）に示したようなエネルギー関数Ｅ（ｘ）である。
【００３２】
以下では、複数の状態変数の数（状態変数の総数）をＮ個とし、Ｎ個の状態変数を状態変数ｘｉ（ｉ＝１〜Ｎ）または状態変数ｘ１〜ｘＮと表記する。また、状態変数ｘ１〜ｘＮの各値の組合せを状態ｘという。なお、以下の説明では、状態保持部１１は、式（２）に示したローカルフィールドｈｉ（ｉ＝１〜Ｎ）についても保持するものとして説明するが、ローカルフィールドｈｉは、エネルギー変化計算部１２が保持していてもよい。
【００３３】
状態保持部１１は、たとえば、レジスタやＳＲＡＭ（Static Random Access Memory）などである。
エネルギー変化計算部１２は、状態変数ｘ１〜ｘＮの値の何れかが変化することに応じて状態遷移が起こる場合、更新インデックス値に基づいて選択された重み値に基づき、エネルギーの変化値（以下エネルギー変化ΔＥｉという）を各状態遷移に対して計算する。
【００３４】
エネルギー変化ΔＥｉは、前述の式（２）のように表せ、エネルギー変化計算部１２は、Δｘｉに応じてローカルフィールドｈｉに符号（＋１または−１）を乗じることでエネルギー変化ΔＥｉを計算できる。また、状態変数ｘｊが０から１に変化したときのローカルフィールドｈｉの変化分Δｈｉは、＋Ｗｉｊ、状態変数ｘｊが１から０に変化したときの変化分Δｈｉは、−Ｗｉｊと表せる。したがって、ローカルフィールドｈｉは行列演算により毎回計算しなおす必要はなく、状態遷移にともなって変化のあったビットによる変化分だけＷｉｊを加減算すればよい。エネルギー変化計算部１２は、計算（または更新）したローカルフィールドｈｉを状態保持部１１に保持させる。
【００３５】
このようなエネルギー変化計算部１２は、たとえば、重み値Ｗｉｊなどを記憶する記憶部（レジスタやＳＲＡＭなど）、セレクタ、乗算器、加算器などを用いて実現できる。なお、以下の説明ではエネルギー変化計算部１２は、Ｎ個のエネルギー変化ΔＥｉ（ｉ＝１〜Ｎ）（以下エネルギー変化ΔＥ１〜ΔＥＮと表記する場合もある）とともに、インデックス値＝ｉ（１〜Ｎ）についても出力するものとする。
【００３６】
オフセット制御部１３は、エネルギー変化ΔＥ１〜ΔＥＮに基づいて各状態遷移のうちの少なくとも１つが許容されるようにオフセット値Ｅｏｆｆを決定し、エネルギー変化ΔＥ１〜ΔＥＮのそれぞれにオフセット値Ｅｏｆｆを付加した複数の評価値を出力する。また、オフセット制御部１３は、所定試行回数毎にオフセット値Ｅｏｆｆを出力する。オフセット制御部１３の例やオフセット値Ｅｏｆｆの決定方法の例については後述する。
【００３７】
比較部１４は、温度値Ｔと乱数値とに基づいて決まる閾値と、複数の評価値のそれぞれとを比較した結果に基づき、各状態遷移（の候補）について状態遷移を許容するか否かを示す複数のフラグ値Ｆ１，Ｆ２，…，ＦＮを出力する。比較部１４は、このような処理によって、式（３）に示した許容確率Ａ（ΔＥ）で状態遷移を許容させることができる。比較部１４の例や、許容確率Ａ（ΔＥ）で状態遷移を許容させることができる理由については後述する。なお、以下の説明では、サンプリング動作は一定の温度で行われる（温度値Ｔが固定値となる）ものとするが、温度は変更されるようにしてもよい。
【００３８】
フラグ計数部１５は、フラグ値Ｆ１〜ＦＮのうち、状態遷移を許容することを示すフラグ値の数を計数し、所定試行回数毎に計数したフラグ値の数（計数値Ｎｆ）を出力する。フラグ計数部１５の例については後述する。
【００３９】
選択部１６はフラグ値Ｆ１〜ＦＮの何れかを選択するとともに、選択した何れかのフラグ値に対応するインデックス値を更新インデックス値として出力する。選択部１６の例についても後述する。
【００４０】
制御部１７は、サンプリング装置１０の各部を制御する。制御部１７は、たとえば、試行回数が所定試行回数に達したか否かを判定し、所定試行回数に達したときに、そのときの状態を状態保持部１１に出力させ、オフセット値Ｅｏｆｆをオフセット制御部１３に出力させる。さらに、制御部１７は、フラグ計数部１５にそのときの計数値Ｎｆを計数させる。また、制御部１７は、比較部１４に温度値Ｔを供給する。
【００４１】
制御部１７は、たとえば、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの特定用途の電子回路にて実現できる。なお、制御部１７は、ＣＰＵ（Central Processing Unit）やＤＳＰ（Digital Signal Processor）などのプロセッサであってもよい。その場合、プロセッサは、図示しないメモリに記憶されたプログラムを実行することで、上記の処理を行う。
【００４２】
試行回数計算部１８は、オフセット制御部１３が出力したオフセット値Ｅｏｆｆと、フラグ計数部１５が出力した計数値Ｎｆに基づき、通常のＭＣＭＣ法において、１つの状態に滞在する試行回数Ｎｔｒｉａｌの期待値＜Ｎｔｒｉａｌ＞（近似値）を計算する。試行回数計算部１８は、たとえば、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路、ＣＰＵやＤＳＰなどのプロセッサにて実現できる。なお、試行回数計算部１８は、サンプリング装置１０の外部に設けられていてもよいし、制御部１７が、オフセット値Ｅｏｆｆと、計数値Ｎｆに基づいて、期待値＜Ｎｔｒｉａｌ＞を計算してもよい。
【００４３】
オフセット値Ｅｏｆｆと計数値Ｎｆを、期待値＜Ｎｔｒｉａｌ＞を求めるための値として用いることができる理由は以下の通りである。
図１には、通常のＭＣＭＣ法における遷移確率Γｉの分布の例が示されている。通常のＭＣＭＣ法では、遷移確率Γｉを全状態変数について積算した値（試行毎に別の状態に遷移する確率（式（４））は、１に満たない。前述のように、式（４）で表される確率の逆数が、通常のＭＣＭＣ法において、１つの状態に滞在する試行回数Ｎｔｒｉａｌの期待値＜Ｎｔｒｉａｌ＞となる。期待値＜Ｎｔｒｉａｌ＞は、現在の状態に留まっている時間（滞在時間）の期待値に相当する。
【００４４】
一方、サンプリング装置１０において、各状態に留まる試行回数を１回とした場合、各状態変数が反転する対象として選ばれる選択確率の分布は、たとえば、以下のようになる。
【００４５】
図２は、各状態に留まる試行回数を１回とした場合の選択確率の分布例を示す図である。
ｐ１〜ｐＮは、状態変数ｘｉ（ｉ＝１〜Ｎ）の選択確率である。各状態に留まる試行回数を１回とした場合、各試行において状態変数ｘｉ（ｉ＝１〜Ｎ）の何れかの変化が生じるため、ｐ１〜ｐＮの積算結果は１となる。選択確率ｐｉは、以下の式（６）で表せる。
【００４６】
【数６】
2020135727.tif 000008
【００４７】
式（６）の右辺は、式（５）と同じになる。
サンプリング装置１０は、各状態遷移のうちの少なくとも１つが許容されるようなオフセット値Ｅｏｆｆを用いることで、上記のように各状態に留まる試行回数が１回となるようにする。オフセット値Ｅｏｆｆを適用したときの許容確率をＡ（ΔＥｉ−Ｅｏｆｆ）とすると、Ａ（ΔＥｉ−Ｅｏｆｆ）は、前述の式（３）により（メトロポリス法を用いた場合）、以下の式（７）のように表せる。
【００４８】
【数７】
2020135727.tif 000009
【００４９】
逆温度βとオフセット値Ｅｏｆｆが正の場合、許容確率Ａ（ΔＥｉ−Ｅｏｆｆ）は、本来の許容確率Ａ（ΔＥｉ）に対して、１以上の係数ｅｘｐ（β・Ｅｏｆｆ）を乗じた確率となる。このため、許容確率Ａ（ΔＥｉ−Ｅｏｆｆ）が用いられる場合、本来の許容確率Ａ（ΔＥｉ）を用いた場合に対して、各状態遷移の許容確率の相互の比率は変わらない。
【００５０】
期待値＜Ｎｔｒｉａｌ＞は、式（７）を用いて、以下の式（８）のように表せる。
【００５１】
【数８】
2020135727.tif 000010
【００５２】
すなわち、状態変数ｘ１〜ｘＮについてのＡ（ΔＥｉ−Ｅｏｆｆ）を積算した値＜Ｎｆ＞は、状態遷移を許容することを示すフラグ値の数である計数値Ｎｆにより近似できる。
式（８）において、逆温度βは一定値（温度値Ｔが固定値であるため）であり、状態変数の総数＝Ｎも固定値である。そのため、サンプリングタイミングにおけるオフセット値Ｅｏｆｆと計数値Ｎｆが得られれば、そのサンプリングタイミングにおける期待値＜Ｎｔｒｉａｌ＞を近似的に求めることができる。
【００５３】
サンプリング装置１０は、オフセット値Ｅｏｆｆを用いて各状態に留まる試行回数が１回となるようにして、所定試行回数毎、状態と、オフセット値Ｅｏｆｆと、計数値Ｎｆとを出力（サンプリング）する。これにより、上記のようにオフセット値Ｅｏｆｆと、計数値Ｎｆとから得られる期待値＜Ｎｔｒｉａｌ＞による重み付けにより、得られた状態が通常のＭＣＭＣ法ではどのくらいの試行回数、その状態に留まるのかがわかり、通常のＭＣＭＣ法の確率過程を再現できる。
【００５４】
図１には、通常のＭＣＭＣ法を用いた場合と、サンプリング装置１０を用いた場合のエネルギーＥと試行回数との関係の一例を示すグラフが示されている。図１に示すように、たとえば、あるエネルギーＥａをもつ状態ｘの滞在時間は、サンプリング装置１０を用いた場合は１試行回数分であるのに対し、通常のＭＣＭＣ法を用いた場合、期待値＜Ｎｔｒｉａｌ＞分である。期待値＜Ｎｔｒｉａｌ＞がわかれば、上記のように通常のＭＣＭＣ法の確率過程を再現できる。
【００５５】
つまり、通常のＭＣＭＣ法にしたがったサンプリングを再現できるため、従来の最適化装置をサンプリング装置として用いた場合に比べて高精度なサンプリングが可能になり、ボルツマン分布にしたがったサンプルが得られる。
【００５６】
機械学習の理論は、サンプルがボルツマン分布になっていることを根拠にしているものがほとんどであり、上記のようなサンプリング装置１０を用いることで、機械学習に適したサンプルを得ることができる。
【００５７】
また、各試行で状態遷移が生じることと、１回のサンプリングで実効的に複数回サンプリングするのと同じ精度が得られることにより、従来の最適化装置をサンプリング装置として用いた場合に比べて計算処理の高速性が損なわれることもない。
【００５８】
以下、サンプリング装置１０の各部の例を説明する。
（状態保持部１１及びエネルギー変化計算部１２の例）
図３は、状態保持部及びエネルギー変化計算部の一例を示す図である。
【００５９】
状態保持部１１は、状態変数ｘ１〜ｘＮの値とローカルフィールドｈ１〜ｈＮとを保持するレジスタ１１ａを有する。レジスタ１１ａは、クロック信号ｃｌｋに同期して、エネルギー変化計算部１２が更新した、状態変数ｘ１〜ｘＮの何れかの値とローカルフィールドｈ１〜ｈＮとを取り込む。クロック信号ｃｌｋは、図示しないクロック信号生成回路から供給される。状態保持部１１は、たとえば、図１に示した制御部１７の制御のもと、所定試行回数毎に複数の状態変数の値を、サンプリング装置１０の外部に出力する。
【００６０】
エネルギー変化計算部１２は、ΔＥ計算回路１２ａ１，１２ａ２，…，１２ａＮ、更新回路１２ｂを有する。
ΔＥ計算回路１２ａ１〜１２ａＮのそれぞれは、状態変数ｘ１〜ｘＮの値の何れかとローカルフィールドｈ１〜ｈＮの何れかを状態保持部１１から読み出し、エネルギー変化ΔＥ１〜ＥＮの何れかを計算する。たとえば、ΔＥ計算回路１２ａ１は、状態変数ｘ１の値とローカルフィールドｈ１とを状態保持部１１から読み出し、状態変数ｘ１の値が変化することによるエネルギー変化ΔＥ１を計算する。ΔＥ計算回路１２ａＮは、状態変数ｘＮの値とローカルフィールドｈＮとを状態保持部１１から読み出し、状態変数ｘＮの値が変化することによるエネルギー変化ΔＥＮを計算する。
【００６１】
エネルギー変化ΔＥｉ（ｉ＝１〜Ｎ）は、式（２）のように表せるため、ΔＥ計算回路１２ａ１〜１２ａＮのそれぞれは、たとえば、状態変数ｘ１〜ｘＮの値に応じてローカルフィールドｈ１〜ｈＮに符号（＋１または−１）を乗じる乗算器を用いて実現できる。
【００６２】
なお、ΔＥ計算回路１２ａ１〜１２ａＮのそれぞれは、インデックス値＝ｉ（１〜Ｎ）の何れかを出力する。たとえば、ΔＥ計算回路１２ａ１は、インデックス値として１を出力し、ΔＥ計算回路１２ａＮは、インデックス値としてＮを出力する。ΔＥ計算回路１２ａ１〜１２ａＮは、インデックス値を保持する保持部（たとえば、レジスタ）を有する。
【００６３】
更新回路１２ｂは、選択部１６より更新インデックス値を受けると、更新インデックス値に対応した状態変数の値を変化させるとともに、その変化に基づいて、ローカルフィールドｈ１〜ｈＮを更新する。更新後の各値は、状態保持部１１に書き込まれる。
【００６４】
たとえば、更新回路１２ｂは、更新インデックス値がｊの場合、状態変数ｘｊを１−ｘｊに変える。そして、更新回路１２ｂは、ローカルフィールドｈｉ（ｉ＝１〜Ｎ）を、ｈｉ＋Ｗｉｊ・Δｘｊに更新する。
【００６５】
なお、更新回路１２ｂは、上記のような更新処理を行うために、現在の状態変数ｘｉ（ｉ＝１〜Ｎ）の値とローカルフィールドｈｉを、状態保持部１１から取得する。重み値Ｗｉｊは、図示しない記憶部に記憶されており、更新回路１２ｂは、その記憶部から重み値Ｗｉｊを取得する。
【００６６】
このような更新回路１２ｂは、重み値Ｗｉｊなどを記憶する記憶部（レジスタやＳＲＡＭなど）、セレクタ、乗算器、加算器などを用いて実現できる。
（オフセット制御部１３及び比較部１４の例）
図４は、オフセット制御部及び比較部の一例を示す図である。なお、図４では、前述のエネルギー変化ΔＥ１〜ΔＥＮがまとめて｛ΔＥｉ｝と表記され、前述のフラグ値Ｆ１〜ＦＮがまとめて｛Ｆｉ｝と表記されている。
【００６７】
オフセット制御部１３は、最小値検出回路１３ａ、加算回路１３ｂ、選択信号生成回路１３ｃ、セレクタ１３ｄを有する。
最小値検出回路１３ａは、エネルギー変化ΔＥ１〜ΔＥＮのうちの最小値ΔＥｍｉｎを検出する。
【００６８】
加算回路１３ｂは、エネルギー変化ΔＥ１〜ΔＥＮから、最小値ΔＥｍｉｎをオフセット値として差し引く（エネルギー変化ΔＥ１〜ΔＥＮに−ΔＥｍｉｎを加算する、ということもできる）演算回路である。加算回路１３ｂの代りに減算回路を用いることもできる。
【００６９】
選択信号生成回路１３ｃは、選択信号ｓｉｇｎを生成する。選択信号生成回路１３ｃは、最小値ΔＥｍｉｎが０より小さい場合、セレクタ１３ｄにエネルギー変化ΔＥ１〜ΔＥＮを選択させる選択信号ｓｉｇｎ（以下では、０とする）を生成する。選択信号生成回路１３ｃは、最小値ΔＥｍｉｎが０以上の場合、セレクタ１３ｄに加算回路１３ｂが算出した値を選択させる選択信号ｓｉｇｎ（以下では、１とする）を生成する。
【００７０】
セレクタ１３ｄは、選択信号ｓｉｇｎが０の場合、エネルギー変化ΔＥ１〜ΔＥＮ（｛ΔＥｉ｝）を選択し、選択信号ｓｉｇｎが１の場合、加算回路１３ｂが算出した値、すなわち｛ΔＥｉ｝−ΔＥｍｉｎを選択して、評価値として出力する。
【００７１】
図４に示すようなオフセット制御部１３を用いた場合、オフセット値Ｅｏｆｆは、Ｅｏｆｆ＝ｍａｘ［０，ｍｉｎ（ΔＥｉ）］と表せる。つまり、エネルギー変化ΔＥ１〜ΔＥＮのうちの最小値ΔＥｍｉｎが０以上の場合、Ｅｏｆｆ＝ΔＥｍｉｎ、それ以外の場合、Ｅｏｆｆ＝０である。
【００７２】
このようなオフセット値Ｅｏｆｆを用いることで、状態変数ｘ１〜ｘＮのうち、少なくとも１つは状態遷移が許容される。その理由を以下に示す。
図５は、オフセット値の適用例を示す図である。図５において、横軸は状態遷移の候補となる状態変数ｘｉを表し、縦軸はエネルギー変化ΔＥを表す。
【００７３】
波形２０ａは、状態変数ｘ１〜ｘＮのそれぞれが変化することによるエネルギー変化｛ΔＥｉ｝の例を示している。図５の例では、エネルギー変化｛ΔＥｉ｝のうち、状態変数ｘｍｉｎの値が変化することによるエネルギー変化ΔＥが最小値ΔＥｍｉｎであり、０より大きい値である。波形２０ｂは、波形２０ａで示されるエネルギー変化｛ΔＥｉ｝から最小値ΔＥｍｉｎを引いた場合のエネルギー変化ΔＥを示す。
【００７４】
波形２０ｂにおけるエネルギー変化ΔＥの最小値は０となるため、状態変数ｘｍｉｎの値が変化する状態遷移の許容確率Ａ（ΔＥ）が１となる。つまり、状態変数ｘ１〜ｘＮのうち、少なくとも状態変数ｘｍｉｎは状態遷移が許容される。
【００７５】
このように、エネルギー変化ΔＥ１〜ΔＥＮのうちの最小値ΔＥｍｉｎが０より大きい場合でも、Ｅｏｆｆ＝ΔＥｍｉｎをエネルギー変化｛ΔＥｉ｝から差し引くことで、状態変数ｘ１〜ｘＮのうち、少なくとも１つの状態遷移が許容される。
【００７６】
図４の説明に戻る。
比較部１４は、乱数発生回路１４ｂ１、選択法則適用部１４ｂ２、乗算回路１４ｂ３、符号反転回路１４ｂ４、比較回路１４ｂ５を有する。
【００７７】
乱数発生回路１４ｂ１は、０より大きく、１以下の一様乱数である乱数値ｒを発生する。乱数発生回路１４ｂ１は、たとえば、メルセンヌツイスタ、ＬＦＳＲ（Linear Feedback Shift Register）などを用いて実現できる。
【００７８】
選択法則適用部１４ｂ２は、選択法則（前述のメトロポリス法またはギブス法）に基づいた値を出力する。
前述の式（３）の許容確率Ａ（ΔＥ）で、エネルギー変化ΔＥを引き起こす状態遷移を許容することを示すフラグ（＝１）を出力する回路は、許容確率Ａ（ΔＥ）と、乱数値ｒとの比較結果に基づいた値を出力する比較器によって実現できる。
【００７９】
ただ、次のような変形を行っても同じ機能が実現できる。２つの数に同じ単調増加関数を作用させても大小関係は変化しない。したがって比較器の２つの入力に同じ単調増加関数を作用させても比較器の出力は変わらない。たとえば、許容確率Ａ（ΔＥ）を、ｆ（−β・ΔＥ）と表記する。ｆ（−β・ΔＥ）に作用させる単調増加関数としてｆ（−β・ΔＥ）の逆関数ｆ−１（−β・ΔＥ）、乱数値ｒに作用させる単調増加関数としてｆ−１（−β・ΔＥ）の−β・ΔＥをｒとしたｆ−１（ｒ）を用いることができる。その場合、上記の比較器と同様の機能を有する回路は、−β・ΔＥがｆ−１（ｒ）より大きいとき１を出力する回路でよいことがわかる。さらにβ＝１／Ｔで、温度値Ｔが正であることから、その回路は、−ΔＥがＴ・ｆ−１（ｒ）より大きいとき１を出力する回路でよい。
【００８０】
選択法則適用部１４ｂ２は、入力される乱数値ｒを上記のｆ−１（ｒ）の値に変換する変換テーブルを用いて、ｆ−１（ｒ）の値を出力する。メトロポリス法が適用される場合、ｆ−１（ｒ）は、ｌｏｇ（ｒ）である。変換テーブルは、たとえば、ＲＡＭ（Random Access Memory）、フラッシュメモリなどのメモリに記憶されている。
【００８１】
乗算回路１４ｂ３は、制御部１７から供給される温度値Ｔと、ｆ−１（ｒ）との積（Ｔ・ｆ−１（ｒ））を閾値として出力する。Ｔ・ｆ−１（ｒ）は、熱励起エネルギーに相当する。
【００８２】
符号反転回路１４ｂ４は、評価値である｛ΔＥｉ｝または｛Ｅｉ｝−ΔＥｍｉｎの符号を反転し、−｛ΔＥｉ｝または−｛ΔＥｉ｝＋ΔＥｍｉｎを出力する。
比較回路１４ｂ５は、−｛ΔＥｉ｝または−｛ΔＥｉ｝＋ΔＥｍｉｎと、Ｔ・ｆ−１（ｒ）とを比較する。そして、比較回路１４ｂ５は、−｛ΔＥｉ｝または−｛ΔＥｉ｝＋ΔＥｍｉｎがＴ・ｆ−１（ｒ）より大きい場合、フラグ値｛Ｆｉ｝として１を出力し、−｛ΔＥｉ｝または−｛ΔＥｉ｝＋ΔＥｍｉｎがＴ・ｆ−１（ｒ）以下の場合、フラグ値｛Ｆｉ｝として０を出力する。
【００８３】
（フラグ計数部１５の例）
図１に示したフラグ計数部１５の例として、以下に、２種のフラグ計数部１５ａ，１５ｂを示す。
【００８４】
図６は、フラグ計数部の一例を示す図である。
フラグ計数部１５ａは、ツリー状に複数段、配置された複数の加算回路（加算回路１５ａ１，１５ａ２，１５ａ３，１５ａ４，１５ａ５など）を有する。
【００８５】
初段の各加算回路は、２つのフラグ値を加算した加算結果を出力する。たとえば、加算回路１５ａ１は、フラグ値Ｆ１，Ｆ２を加算した加算結果を出力し、加算回路１５ａ２は、フラグ値Ｆ３，Ｆ４を加算した加算結果を出力する。２段目以降の加算回路は、前段の２つの加算回路の出力値を加算した加算結果を出力する。たとえば、２段目の加算回路１５ａ４は、初段の加算回路１５ａ１，１５ａ２の出力値を加算した加算結果を出力する。そして、最後段の加算回路１５ａ５が出力する加算結果が、計数値Ｎｆとなる。
【００８６】
初段の各加算回路は、２つの１ビットのフラグ値を加算して、２ビットの加算結果を出力し、２段目の各加算回路は、２つの加算回路のそれぞれの２ビットの出力値を加算して、３ビットの加算結果を出力する。Ｎ＝１０２４の場合、最後段（１０段目）の加算回路１５ａ５は、前段の２つの加算回路のそれぞれの１０ビットの出力値を加算して、１１ビットの計数値Ｎｆを出力する。
【００８７】
図７は、フラグ計数部の他の例を示す図である。
フラグ計数部１５ｂは、マイクロコントローラ１５ｂ１とプログラムメモリ１５ｂ２とを有する。
【００８８】
マイクロコントローラ１５ｂ１は、プログラムメモリ１５ｂ２に記憶されたプログラムを読み出して実行することで、フラグ値Ｆ１〜ＦＮのうち値が１である数を計数し、その計数結果である計数値Ｎｆを出力する。
【００８９】
サンプリング装置１０によるサンプリングは、通常１０００クロックサイクル以上の間隔（サンプリング間隔）で行われるため、比較的処理速度が遅いマイクロコントローラ１５ｂ１及び、比較的小規模のプログラムメモリ１５ｂ２を用いることができる。
【００９０】
（選択部１６の例）
図８は、選択部の一例を示す図である。
選択部１６は、乱数発生回路１６ａ、加算回路１６ｂ、選択回路部１６ｃを有する。
【００９１】
乱数発生回路１６ａは、たとえば、０＜ｒｉ＜１の一様乱数である乱数値ｒｉを、インデックス値＝ｉ（１〜Ｎ）のそれぞれについて異なるシードを用いて発生する。乱数発生回路１６ａは、たとえば、メルセンヌツイスタ、ＬＦＳＲなどを用いて実現できる。
【００９２】
加算回路１６ｂは、フラグ値Ｆ１〜ＦＮのそれぞれに対して、複数の乱数値ｒｉ（ｉ＝１〜Ｎ）の何れかを加算した加算結果を、フラグ値Ｆ１〜ＦＮに対応するインデックス値とともに出力する。
【００９３】
選択回路部１６ｃは、ツリー状に複数段、配置された複数の選択回路（選択回路１６ｃ１，１６ｃ２，１６ｃ３，１６ｃ４，１６ｃ５など）を有する。
初段の各選択回路は、２つのフラグ値についての加算結果のうち、大きい方の加算結果とそれに対応したインデックス値を選択して出力する。たとえば、選択回路１６ｃ１は、Ｆ１＋ｒ１がＦ２＋ｒ２より大きい場合、Ｆ１＋ｒ１と、インデックス値＝１を出力し、Ｆ１＋ｒ１がＦ２＋ｒ２より小さい場合、Ｆ２＋ｒ２と、インデックス値＝２を出力する。
【００９４】
２段目以降の選択回路は、前段の２つの選択回路が選択した加算結果のうち、大きい方の加算結果とそれに対応したインデックス値を選択して出力する。たとえば、２段目の選択回路１６ｃ４は、初段の選択回路１６ｃ１，１６ｃ２が選択した加算結果のうち、大きい方の加算結果とそれに対応したインデックス値を選択して出力する。最後段の選択回路１６ｃ５は、前段の２つの選択回路が選択した加算結果のうち、大きい方の加算結果に対応したインデックス値＝ｊを更新インデックス値として出力する。
【００９５】
このように、選択部１６は、加算結果の大小関係に基づいて、フラグ値Ｆ１〜ＦＮの何れかを選択することで、フラグ値Ｆ１〜ＦＮのうち、値が１のものを、ランダムに等確率で選択することができる。
【００９６】
（動作例）
以下、サンプリング装置１０の動作例を説明する。
なお、以下の例ではサンプリング動作は、温度値Ｔが一定の値で行われるものとする。
【００９７】
図９は、サンプリングタイミングの例を示す図である。
図９において、期間Ｔｕは、１回の試行（状態更新処理）にかかる期間を示す。サンプリング装置１０は、所定の試行回数（バーンイン期間）の試行が行われたのちに、サンプリングを開始する。バーンイン期間は、状態が平衡状態に達するまでの期間であり、バーンイン期間に対応する試行回数が予め設定される。
【００９８】
サンプリング期間では、サンプリング装置１０は、所定の試行回数（サンプリング間隔Ｔｓ）毎に、サンプリング（たとえば、状態ｘ、オフセット値Ｅｏｆｆ、計数値Ｎｆなどの出力）を行う。サンプリング間隔Ｔｓに対応する試行回数は予め設定される。
【００９９】
図１０は、サンプリング動作の一例の流れを示すフローチャートである。
まず、初期化やパラメータ設定が行われる（ステップＳ１）。ステップＳ１の処理において、制御部１７は、たとえば、重み値、バイアス係数、状態変数ｘｉ（ｉ＝１〜Ｎ）の初期値に基づいて、ローカルフィールドｈｉ（ｉ＝１〜Ｎ）の初期値の計算を行う。重み値、バイアス係数、状態変数ｘｉ（ｉ＝１〜Ｎ）の初期値は、計算対象の組合せ最適化問題を変換したイジングモデルの情報として、たとえば、図示しない記憶部に予め記憶されている。状態変数ｘｉの初期値やローカルフィールドｈｉの初期値は、状態保持部１１に保持される。
【０１００】
また、バーンイン期間に対応する試行回数ＮＢ、サンプリング間隔Ｔｓに対応する試行回数ＮＳＩ、全サンプルが得られる試行回数ＮＳＮについても、たとえば、図示しない記憶部に予め記憶されている。そして、ステップＳ１の処理において、制御部１７にそれらのパラメータが設定される。
【０１０１】
その後、エネルギー変化計算部１２は、状態変数ｘｉの値とローカルフィールドｈｉを、状態保持部１１から読み出す（ステップＳ２）。
そして、エネルギー変化計算部１２は、前述のように、エネルギー変化ΔＥｉ（ｉ＝１〜Ｎ）を計算する（ステップＳ３）。
【０１０２】
オフセット制御部１３は、エネルギー変化計算部１２から供給されるエネルギー変化ΔＥｉに対して、前述のようなオフセット値Ｅｏｆｆを付加する（ステップＳ４）。
そして、前述の比較部１４と選択部１６の処理により、状態遷移を許容する状態変数を識別する更新インデックス値が出力される（ステップＳ５）。
【０１０３】
エネルギー変化計算部１２は、更新インデックス値により識別される状態変数ｘｉの値を更新する。また、エネルギー変化計算部１２は、その更新に基づいて、ローカルフィールドｈｉを、前述の処理により更新する（ステップＳ６）。
【０１０４】
更新後の状態変数ｘｉとローカルフィールドｈｉの各値は、状態保持部１１に書き込まれる（ステップＳ７）。
制御部１７は、ステップＳ２〜Ｓ７の処理（試行）の回数が試行回数ＮＢに達したか否かを判定する（ステップＳ８）。ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＢに達していない場合、ステップＳ２からの処理が繰り返される。
【０１０５】
ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＢに達した場合、制御部１７は、ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＳＩに達したか否かを判定する（ステップＳ９）。ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＳＩに達していない場合、ステップＳ２からの処理が繰り返される。
【０１０６】
ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＳＩに達した場合、制御部１７は、フラグ計数部１５に、計数値Ｎｆを計数させる（ステップＳ１０）。そして、制御部１７は、フラグ計数部１５に計数値Ｎｆ、状態保持部１１に状態ｘ、オフセット制御部１３にオフセット値Ｅｏｆｆをそれぞれ出力させる（ステップＳ１１）。出力されたオフセット値Ｅｏｆｆと、計数値Ｎｆに基づいて、試行回数計算部１８は、式（８）に示される期待値＜Ｎｔｒｉａｌ＞を計算する。
【０１０７】
制御部１７は、ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＳＮに達したか否かを判定する（ステップＳ１２）。ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＳＮに達していない場合、ステップＳ２からの処理が繰り返される。ステップＳ２〜Ｓ７の処理の回数が試行回数ＮＳＮに達した場合、制御部１７は、サンプリング動作を終了させる。
【０１０８】
ところで、サンプリング装置１０は、最適化装置としても使用できる。最適化装置として使用する場合の動作例を以下に示す。
図１１は、最適化処理動作の一例の流れを示すフローチャートである。
【０１０９】
まず、初期化やパラメータ設定が行われる（ステップＳ２０）。ステップＳ２０の処理において、制御部１７は、サンプリング動作時と同様に、ローカルフィールドｈｉ（ｉ＝１〜Ｎ）の初期値の計算などを行う。
【０１１０】
また、ステップＳ２０の処理において、たとえば、図示しない記憶部に予め記憶された温度変更スケジュールにしたがって、制御部１７は、温度値Ｔの初期値を比較部１４に設定する。また、たとえば、図示しない記憶部には予め、温度値Ｔの更新を行う試行回数Ｎ１、上限の試行回数Ｎ２が記憶されており、ステップＳ２０の処理において、制御部１７にそれらのパラメータが設定される。
【０１１１】
その後のステップＳ２１〜Ｓ２６の処理は、図１０に示したサンプリング動作時におけるステップＳ２〜Ｓ７の処理とほぼ同様である。ただし、たとえば、制御部１７は、式（１）で表せるエネルギーＥの初期値を計算するとともに、更新インデックス値に対応するエネルギー変化（たとえば、更新インデックス値＝ｊの場合、エネルギー変化ΔＥｊ）を用いて、エネルギーＥを更新する。そして、制御部１７は、これまでの最低エネルギーとそのときの状態を図示しない記憶部に保持させる。
【０１１２】
ステップＳ２６の後、制御部１７は、ステップＳ２１〜Ｓ２６の処理（試行）の回数が試行回数Ｎ１に達したか否かを判定する（ステップＳ２７）。ステップＳ２１〜Ｓ２６の処理の回数が試行回数Ｎ１に達していない場合、ステップＳ２１からの処理が繰り返される。
【０１１３】
ステップＳ２１〜Ｓ２６の処理の回数が試行回数Ｎ１に達した場合、制御部１７は、シミュレーテッド・アニーリングを実現するために、温度値Ｔを、温度変更スケジュールにしたがって小さくするように更新する（ステップＳ２８）。
【０１１４】
その後、制御部１７は、ステップＳ２１〜Ｓ２６の処理の回数が試行回数Ｎ２に達したか否かを判定する（ステップＳ２９）。ステップＳ２１〜Ｓ２６の処理の回数が試行回数Ｎ２に達していない場合、ステップＳ２１からの処理が繰り返される。ステップＳ２１〜Ｓ２６の処理の回数が試行回数Ｎ２に達した場合、制御部１７は、その時点で記憶部に保持させている最低エネルギーのときの状態を、組合せ最適化問題の解として出力し（ステップＳ３０）、最適化処理動作を終了させる。
【０１１５】
サンプリング装置１０を最適化装置として用いた場合、前述のようなオフセット値Ｅｏｆｆを用いることで各試行において状態遷移が生じ、探索が高速化する。
（第２の実施の形態）
以下、第２の実施の形態のサンプリング装置を説明する。
【０１１６】
前述のように値が１であるフラグ値の数を数えることは、許容確率Ａ（ΔＥｉ）（ｉ＝１〜Ｎ）を積算したものを、仮数部１ビット、指数部をｅｘｐ（−β・Ｅｏｆｆ）で表すことに相当する。以下に示す第２の実施の形態のサンプリング装置では、２種類のフラグ値から計数値Ｎｆを求めることで、より精度のよい期待値＜Ｎｔｒｉａｌ＞を計算するものである。
【０１１７】
図１２は、第２の実施の形態のサンプリング装置の一例を示す図である。図１２において、図１に示した要素と同じ要素については同じ符号が付されている。なお、図１２では、図１に示した制御部１７については図示が省略されている。また、以下では、図１においてオフセット制御部１３が出力するオフセット値Ｅｏｆｆをオフセット値Ｅｏｆｆ１といい、図１においてフラグ計数部１５が出力する計数値Ｎｆを計数値Ｎｆ１という。
【０１１８】
第２の実施の形態のサンプリング装置３０は、ΔＥ置き換え部３１、オフセット制御部３２、比較部３３、フラグ計数部３４、選択部３５、計数値補正部３６、試行回数計算部３７を有する。
【０１１９】
ΔＥ置き換え部３１は、フラグ値Ｆ１〜ＦＮに基づき、エネルギー変化ΔＥ１〜ΔＥＮのうち、許容されない状態遷移であるとされている複数の状態遷移に対する複数の非許容エネルギー変化値（以下エネルギー変化ΔＥｎａ１〜ΔＥｎａＭという）を出力する。ΔＥ置き換え部３１は、たとえば、エネルギー変化ΔＥ１〜ΔＥＮと、フラグ値Ｆ１〜ＦＮを入力し、値が０のフラグ値のインデックス値に対応するインデックス値をもつエネルギー変化をエネルギー変化ΔＥ１〜ΔＥＮから選択する回路によって実現できる。
【０１２０】
オフセット制御部３２は、エネルギー変化ΔＥｎａ１〜ΔＥｎａＭに基づき、フラグ値Ｆ１〜ＦＮでは許容されない状態遷移とされている複数の状態遷移のうちの少なくとも１つが許容されるようにオフセット値Ｅｏｆｆ１とは別のオフセット値Ｅｏｆｆ２を決定する。そして、オフセット制御部３２は、エネルギー変化ΔＥｎａ１〜ΔＥｎａＭのそれぞれにオフセット値Ｅｏｆｆ２を付加した、オフセット制御部１３が出力する複数の評価値とは別の複数の評価値を出力する。また、オフセット制御部３２は、所定試行回数毎にオフセット値Ｅｏｆｆ２を出力する。オフセット制御部３２は、オフセット制御部１３と同様の回路にて実現できる。
【０１２１】
比較部３３は、温度値Ｔと乱数値とに基づいて決まる閾値と、比較部３３が出力する複数の評価値のそれぞれとを比較した結果に基づき複数の状態遷移のそれぞれについて状態遷移を許容するか否かを示す、フラグ値Ｆ１〜ＦＮとは別の複数のフラグ値を出力する。以下、比較部３３が出力する複数のフラグ値をフラグ値Ｆｎａ１〜ＦｎａＭという。比較部３３は、比較部１４と同様の回路にて実現できる。
【０１２２】
フラグ計数部３４は、比較部３３が出力するフラグ値Ｆｎａ１〜ＦｎａＭのうち、状態遷移を許容することを示すフラグ値（値が１）の数を計数した、計数値Ｎｆ１とは別の計数値Ｎｆ２を、所定試行回数毎に出力する。フラグ計数部３４は、フラグ計数部１５と同様の回路にて実現できる。
【０１２３】
選択部３５は、計数値Ｎｆ１，Ｎｆ２、オフセット値Ｅｏｆｆ１，Ｅｏｆｆ２に基づき、Ｎｆ１：Ｎｆ２ｅｘｐ［−β（Ｅｏｆｆ１−Ｅｏｆｆ２）］の比率で、フラグ値Ｆ１〜ＦＮのうち値が１のグループと、フラグ値Ｆｎａ１〜ＦｎａＭのうち値が１のグループの何れか一方を選択する。そして、選択部３５は、選択したグループから等確率に１つのフラグ値を選択し、選択したフラグ値に対応したインデックス値を更新インデックス値として出力する。選択部３５は、たとえば、論理回路、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路、ＣＰＵやＤＳＰなどのプロセッサ、マイクロコントローラなどで実現できる。
【０１２４】
計数値補正部３６は、計数値Ｎｆ２とオフセット値Ｅｏｆｆ１，Ｅｏｆｆ２に基づいて、計数値Ｎｆ１を補正し、計数値Ｎｆとして出力する。計数値補正部３６は、たとえば、Ｎｆ＝Ｎｆ１＋Ｎｆ２ｅｘｐ［β（Ｅｏｆｆ１−Ｅｏｆｆ２）］を計算して出力する。計数値補正部３６は、たとえば、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路、ＣＰＵやＤＳＰなどのプロセッサにて実現できる。
【０１２５】
試行回数計算部３７は、オフセット値Ｅｏｆｆ１，Ｅｏｆｆ２と、計数値Ｎｆに基づいて、たとえば、前述の式（８）を用いて期待値＜Ｎｔｒｉａｌ＞を計算する。ただし、式（８）のＥｏｆｆとして、（Ｅｏｆｆ１−Ｅｏｆｆ２）が用いられる。試行回数計算部３７は、たとえば、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路、ＣＰＵやＤＳＰなどのプロセッサにて実現できる。
【０１２６】
その他の構成及び動作については、第１の実施の形態のサンプリング装置１０と同じである。
上記のような第２の実施の形態のサンプリング装置３０によれば、たとえば、フラグ値Ｆ１〜ＦＮのうち、１つだけが１であるような（１つの状態遷移だけが許容されるような）極端な場合でも、他の状態遷移についてもある割合で許容されるようになる。これにより、サンプリング装置３０を用いた場合の各状態遷移（状態変数）間の遷移確率の比を、通常のＭＣＭＣ法で得られるような各状態遷移間の遷移確率の比に近づけることができる。そして、上記のように補正された計数値Ｎｆを用いて計算される期待値＜Ｎｔｒｉａｌ＞は、通常のＭＣＭＣ法で得られる１つの状態に滞在する試行数により近いものになり、通常のＭＣＭＣ法により近い確率過程が得られる。
【０１２７】
なお、上記の第２の実施の形態のサンプリング装置３０の動作は、１つのオフセット制御部１３、１つの比較部１４、１つのフラグ計数部１５を用いても実行できる。その場合、以下のような動作によって、２種類のフラグ値（フラグ値Ｆ１〜ＦＮとフラグ値Ｆｎａ１〜ＦｎａＭ）についての計数値Ｎｆ１，Ｎｆ２が得られる。
【０１２８】
図１３は、２種類のフラグ値についての計数値Ｎｆ１，Ｎｆ２を発生させる動作の一例の流れを示すフローチャートである。
以下に示す処理は、たとえば、図１に示した制御部１７（図１２では図示が省略されている）によって制御される。
【０１２９】
まず、制御部１７は、ｎＭＡＸ＝２、ｎ＝１とする（ステップＳ４０）。そして、制御部１７は、オフセット制御部１３に、オフセット値Ｅｏｆｆｎ（ｎ＝１の場合はオフセット値Ｅｏｆｆ１）を発生させる（ステップＳ４１）。オフセット制御部１３は、エネルギー変化ΔＥｉ（ｉ＝１〜Ｎ）のそれぞれからオフセット値Ｅｏｆｆｎを差し引いた評価値を計算する（ステップＳ４２）。図４のようなオフセット制御部１３が用いられる場合、ΔＥｉの最小値ΔＥｍｉｎが０より小さい場合、オフセット値Ｅｏｆｆｎは０であり、最小値ΔＥｍｉｎが０以上の場合、オフセット値Ｅｏｆｆｎは最小値ΔＥｍｉｎである。
【０１３０】
比較部１４は、温度値Ｔと乱数値とに基づいて決まる閾値と、複数の評価値のそれぞれとを比較した結果に基づき、フラグ値（ｎ＝１の場合、フラグ値Ｆ１〜ＦＮ）を発生する（ステップＳ４３）。
【０１３１】
そして、フラグ計数部１５は、フラグ値が１である数を計数し（ステップＳ４４）、計数値Ｎｆｎ（ｎ＝１の場合は計数値Ｎｆ１）を出力する（ステップＳ４５）。
制御部１７は、ｎ≧ｎＭＡＸであるか否かを判定し（ステップＳ４６）、ｎ≧ｎＭＡＸの場合、計数値Ｎｆ１，Ｎｆ２を発生させる処理を終了する。ｎ＜ｎＭＡＸの場合、制御部１７は、ｎ＝ｎ＋１とする（ステップＳ４７）。また、制御部１７は、前述のΔＥ置き換え部３１の機能により、エネルギー変化ΔＥｉを、エネルギー変化ΔＥｎａ１〜ΔＥｎａＭに置き換え（ステップＳ４８）、ステップＳ４１からの処理を繰り返す。
【０１３２】
なお、上記ステップＳ４４，Ｓ４５の処理以外は、各試行において行われる。
（第３の実施の形態）
図１４は、第３の実施の形態のサンプリング装置の一例を示す図である。図１４において、図１に示した要素と同じ要素については図示が省略されている。
【０１３３】
第３の実施の形態のサンプリング装置４０では、図１に示した選択部１６の代りに、選択部４１が設けられている。
選択部４１は、エネルギー変化ΔＥ１〜ΔＥＮのうち、０より小さいものを０に更新させた複数の更新エネルギー変化値を生成する。そして、選択部４１は、複数の更新エネルギー変化値のそれぞれに、複数の閾値の何れかを加算した複数の評価値のうちで、最小となる評価値に対応するインデックス値を更新インデックス値として出力する。複数の閾値は、温度値Ｔと互いに独立な複数の乱数値とに基づいてそれぞれ算出される。
【０１３４】
選択部４１は、更新エネルギー変化値算出回路４１ａ、閾値生成部４１ｂ、加算回路４１ｃ１，４１ｃ２，…，４１ｃＮ、最小値選択回路４１ｄを有する。
更新エネルギー変化値算出回路４１ａは、エネルギー変化ΔＥ１〜ΔＥＮのうち、０以上のものをそのまま出力し、０より小さいものを０に更新することで、Ｎ個の更新エネルギー変化値を生成する。更新エネルギー変化値算出回路４１ａの処理は、ｍａｘ［０，ΔＥｉ］を計算することに相当する。
【０１３５】
閾値生成部４１ｂは、Ｎ個の閾値を生成する。閾値生成部４１ｂの例については後述する。
加算回路４１ｃ１〜４１ｃＮは、更新エネルギー変化値算出回路４１ａが出力したＮ個の更新エネルギー変化値のそれぞれに、複数の閾値の何れかを加算したＮ個の評価値を出力する。
【０１３６】
最小値選択回路４１ｄは、加算回路４１ｃ１〜４１ｃＮが出力するＮ個の評価値のうちで、最小となる評価値に対応するインデックス値を更新インデックス値として出力する。なお、インデックス値は、たとえば、エネルギー変化計算部１２から供給される。最小値選択回路４１ｄは、たとえば、ツリー状に複数段、配置された複数の選択回路を有し、各選択回路は２つの評価値のうち、小さい方を選択して出力する。最後段の選択回路から出力される評価値に対応したインデックス値が更新インデックス値となる。
【０１３７】
図１５は、閾値生成部の一例を示す図である。
閾値生成部４１ｂは、メルセンヌツイスタ５０ａ１，５０ａ２，…，５０ａｎ、変換部５０ｂ１，５０ｂ２，…，５０ｂｎ、レジスタ５０ｃ１１〜５０ｃｎｍ、乗算部５０ｄを有する。
【０１３８】
メルセンヌツイスタ５０ａ１〜５０ａｎは、それぞれ０＜ｒｉ＜１の一様乱数である乱数値ｒｉを、インデックス値＝ｉ（１〜Ｎ）のそれぞれについて異なるシードを用いて発生する。乱数値ｒｉは、たとえば、１６ビットの値で表される。なお、乱数値ｒｉは、１クロックサイクル毎に更新される。
【０１３９】
変換部５０ｂ１〜５０ｂｎは、変換テーブルを用いて、乱数値ｒｉをｌｏｇ（−ｌｏｇｒｉ）に変換する。ｌｏｇ（−ｌｏｇｒｉ）は、たとえば、２７ビットの値で表される。変換テーブルは、たとえば、ＲＡＭ、フラッシュメモリなどのメモリに記憶されている。
【０１４０】
レジスタ５０ｃ１１〜５０ｃｎｍは、図示しないクロック信号に同期して、変換部５０ｂ１〜５０ｂｎが出力した値を遅延して出力する。たとえば、直列に接続されたレジスタ５０ｃ１１，５０ｃ１２，…，５０ｃ１ｍにより、変換部５０ｂ１が出力した値から、ｍ個の異なる乱数値が生成される。また、直列に接続されたレジスタ５０ｃ２１，５０ｃ２２，…，５０ｃ２ｍにより、変換部５０ｂ２が出力した値から、ｍ個の異なる乱数値が生成される。また、直列に接続されたレジスタ５０ｃｎ１，５０ｃｎ２，…，５０ｃｎｍにより、変換部５０ｂｎが出力した値から、ｍ個の異なる乱数値が生成される。
【０１４１】
乗算部５０ｄは、変換部５０ｂ１〜５０ｂｎが出力した値と、レジスタ５０ｃ１１〜５０ｃｎｍのそれぞれが出力した値に、温度値Ｔを乗じた値を閾値として出力する。乗算部５０ｄは、Ｎ個の乗算回路を有する。
【０１４２】
たとえば、Ｎ＝１０２４の場合、メルセンヌツイスタ５０ａ１〜５０ａｎと変換部５０ｂ１〜５０ｂｎは、それぞれ３２個（ｎ＝３２）、レジスタ５０ｃ１１〜５０ｃｎｍの数（ｎ×ｍ）は、３２×３１である。これにより、１０２４個の閾値が生成される。
【０１４３】
選択部１６の代りに上記の選択部４１を用いることで、試行毎に適切な許容確率で状態遷移が発生する。その理由を以下に説明する。
前述した０＜ｒｉ＜１の一様乱数である乱数値ｒｉの累積分布関数Ｆ（ｒ）＝Ｐｒｏｂ（ｒｉ≦ｒ）の値は、ｒ≦０の場合に０、０＜ｒ＜１の場合にｒ、ｒ≧１の場合に１となる。ここで、乱数値ｒｉから発生させる正の乱数値ｙｉを、ｙｉ＝−ｌｏｇ（ｒｉ）／Ａｉ（Ａｉ＞０）とする。乱数値ｙｉがｙ（＞０）より大きくなる確率Ｐｒｏｂ（ｙｉ≧ｙ）は、Ｐｒｏｂ（ｙｉ≧ｙ）＝Ｐｒｏｂ（ｒｉ≦ｅｘｐ（−Ａｉｙ））＝Ｆ（ｅｘｐ（−Ａｉｙ））＝ｅｘｐ（−Ａｉｙ）と表せる。したがって乱数値ｒｉの確率密度関数ｐ（ｙｉ）は、以下の式（９）のように表せる。
【０１４４】
【数９】
2020135727.tif 000011
【０１４５】
あるインデックス値＝ｉに対する乱数値ｙｉが他の全ての乱数値ｙｊ（ｊ≠ｉ）より小さい確率は、以下の式（１０）のように表せる。
【０１４６】
【数１０】
2020135727.tif 000012
【０１４７】
式（１０）の右辺は、式（５）の右辺と同じ形である。
ここで、Ａｉをエネルギー変化ΔＥｉに対する許容確率とし、Ａｉ＝ｍｉｎ［１，ｅｘｐ（−βΔＥｉ）］とすれば、式（１０）が１となるときのＡｉが求めたい許容確率となる。なお、ｌｏｇ（Ａｉ）／β＝−ｍａｘ［０，ΔＥｉ］となる。
【０１４８】
あるインデックス値＝ｉに対する乱数値ｙｉが他の全ての乱数値ｙｊ（ｊ≠ｉ）より小さいことはｌｏｇ（ｙｉ）／βがｌｏｇ（ｙｊ）／β（ｊ≠ｉ）より小さいことと同じである。ここで、ｌｏｇ（ｙｉ）／β＝ｌｏｇ（−ｌｏｇ（ｒｉ））／β＋ｍａｘ［０，ΔＥｉ］、β＝１／Ｔであるから、ｍａｘ［０，ΔＥｉ］にＴ・ｌｏｇ（−ｌｏｇ（ｒｉ））を加算した値の中で最小のものを選べば、試行毎に適切な許容確率で状態遷移が発生する。
【０１４９】
その他の構成及び動作については、第１の実施の形態のサンプリング装置１０と同じである。
図１６は、温度とエネルギーとの関係についてのシミュレーション結果の例を示す図である。
【０１５０】
また、図１７は、試行回数と最低エネルギー状態が得られる確率との関係についてのシミュレーション結果の例を示す図である。図１６、図１７で扱われている問題は、１６都市の巡回セールスマン問題である。
【０１５１】
図１６では、横軸は温度値Ｔ、縦軸は式（１）に示したエネルギーＥを表している。特性６０ａは、通常のＭＣＭＣ法を用いた場合の温度とエネルギーとの関係を示し、特性６０ｂは、第３の実施の形態のサンプリング装置４０を用いた場合の温度とエネルギーとの関係を示している。また、特性６０ｃは、デジタル回路を用いた従来の最適化装置を用いた場合の温度とエネルギーとの関係を示している。
【０１５２】
図１６に示すように、第３の実施の形態のサンプリング装置４０は、通常のＭＣＭＣ法を用いた場合とほぼ同様の特性が得られている。
図１７では、横軸は試行回数、縦軸は最低エネルギー状態が得られる確率を表している。特性６１ａは、通常のＭＣＭＣ法を用いた場合の試行回数と上記確率との関係を示し、特性６１ｂは、第３の実施の形態のサンプリング装置４０を用いた場合の試行回数と上記確率との関係を示している。また、特性６１ｃは、デジタル回路を用いた従来の最適化装置を用いた場合の試行回数と上記確率との関係を示している。
【０１５３】
図１７に示すように、第３の実施の形態のサンプリング装置４０は、従来の最適化装置や通常のＭＣＭＣ法を用いた場合よりも、高速に（早い試行回数で）最低エネルギー状態が得られることがわかる。
【０１５４】
以上のように、第３の実施の形態のサンプリング装置４０は、理想的に１回の試行で状態遷移が行われ、理論的にも通常のＭＣＭＣ法と同等のサンプル分布が得られるため、高速で高精度のサンプリングが可能となる。
【０１５５】
（変形例）
以下に示すサンプリング装置は、レプリカ交換法を利用したものである。
レプリカ交換法は複数の温度を用いたＭＣＭＣ法を同時に行い、ある試行回数毎に、それぞれ状態のエネルギーを比較し、適切な確率で２つの温度に対する状態（または温度）を交換するという操作を行う方法である。以下の例では、温度を交換するものとして説明するが、状態を交換しても同じことである。
【０１５６】
図１８は、レプリカ交換法を利用したサンプリング装置の例を示す図である。
図１８に示すサンプリング装置７０は、サンプリング処理部７１とレプリカ交換コントローラ７２とを有する。
【０１５７】
サンプリング処理部７１は、図１に示したサンプリング装置１０の各部を有しており、パイプライン処理によって、異なる温度値Ｔ（図１８の例では逆温度β１，β２，…，βＭ）が設定されたレプリカ７１ａ１，７１ａ２，…，７１ａＭのそれぞれについて、前述の試行を行う。
【０１５８】
状態ｘやエネルギーＥは、レプリカ７１ａ１〜７１ａＭのそれぞれについて保持され、試行が行われるたびに更新される。
レプリカ交換コントローラ７２は、所定試行回数毎にレプリカ７１ａ１〜７１ａＭのエネルギーＥを比較して、その比較結果に基づいてレプリカ７１ａ１〜７１ａＭのうちの２つにおいて設定されている逆温度（または温度）を、所定の交換確率ｐｉｊで交換する。交換確率ｐｉｊは、ｐｉｊ＝ｆ（（βｉ−βｊ）（Ｅｉ−Ｅｊ））と表せる。βｉはｉ番目のレプリカに設定された逆温度、βｊはｊ番目のレプリカに設定された逆温度、Ｅｉはｉ番目のレプリカにおけるエネルギー、Ｅｊはｊ番目のレプリカにおけるエネルギーである。また、関数ｆは式（３）のβの代りにβｉ−βｊが用いられ、式（３）のΔＥの代りにＥｉ−Ｅｊが用いられる。
【０１５９】
レプリカ交換コントローラ７２は、たとえば、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路、または、ＣＰＵやＤＳＰなどのプロセッサなどである。
なお、上記のレプリカ交換コントローラ７２の処理を、図１に示した制御部１７が行ってもよい。
【０１６０】
サンプリング動作は、所定の試行回数毎に行われ、たとえば、温度値Ｔが最低のレプリカについての状態ｘと計数値Ｎｆとオフセット値Ｅｏｆｆが出力される。
図１９は、レプリカ交換法を利用したサンプリング装置の動作例を示すフローチャートである。なお、図１０に示した第１の実施の形態のサンプリング装置１０の動作と同じ処理については、図示及び説明を省略する。
【０１６１】
レプリカ交換法を利用する場合、初期化（ステップＳ５０）の際に、レプリカ交換コントローラ７２に、最高逆温度βｍａｘ、最低逆温度βｍｉｎ、交換間隔を示す試行回数Ｎｉ＿ｔｏｔ、総サンプリング回数Ｎｓ＿ｔｏｔが設定される。レプリカ交換コントローラ７２は、これらのパラメータを、たとえば、図示しない記憶部から読み出してもよい。また、初期化の際に、レプリカ交換コントローラ７２は、最高逆温度βｍａｘ、最低逆温度βｍｉｎに基づいて、レプリカ７１ａ１〜７１ａＭに逆温度β１〜βＭの初期値を設定する。
【０１６２】
次に、たとえば、レプリカ交換コントローラ７２は、図示しないカウンタにサンプリング回数Ｎｓとして０を設定し（ステップＳ５１）、別の図示しないカウンタに試行回数Ｎｉ＝０を設定する（ステップＳ５２）。その後、試行が行われ各レプリカの状態ｘとエネルギーＥが更新される（ステップＳ５３）。レプリカ交換コントローラ７２は、試行回数Ｎｉを１増加させ（ステップＳ５４）、Ｎｉ≧Ｎｉ＿ｔｏｔであるか否かを判定する（ステップＳ５５）。Ｎｉ≧Ｎｉ＿ｔｏｔではない場合には、ステップＳ５３からの処理が繰り返される。
【０１６３】
Ｎｉ≧Ｎｉ＿ｔｏｔである場合、レプリカ交換コントローラ７２は、前述の交換確率を用いてレプリカ間の温度（逆温度β）を交換するとともに、サンプリング回数ＮＳを１増加させる（ステップＳ５６）。そして、レプリカ交換コントローラ７２は、温度値Ｔが最低（逆温度βが最高）のレプリカについての計数値Ｎｆと状態ｘとオフセット値Ｅｏｆｆを出力させる（ステップＳ５７）。また、出力された計数値Ｎｆとオフセット値Ｅｏｆｆとに基づいて、図１に示した試行回数計算部１８は、期待値＜Ｎｔｒｉａｌ＞を計算する。
【０１６４】
その後、レプリカ交換コントローラ７２は、Ｎｓ≧Ｎｓ＿ｔｏｔであるか否かを判定する（ステップＳ５８）。Ｎｓ≧Ｎｓ＿ｔｏｔではない場合には、ステップＳ５２からの処理が繰り返される。
【０１６５】
Ｎｓ≧Ｎｓ＿ｔｏｔである場合、レプリカ交換コントローラ７２は、サンプリング動作を終了させる。
このようなレプリカ交換法を利用したサンプリング装置７０は、局所解からの脱出が早く、温度値Ｔの最適化が容易であるという利点がある。
【０１６６】
次に、サンプリング装置の他の変形例を示す。
図２０は、サンプリング装置の他の変形例を示す図である。
図２０に示すサンプリング装置８０は、サンプリング処理部８１と期待値計算部８２を有する。
【０１６７】
サンプリング処理部８１は、図１に示したサンプリング装置１０の各部を有しており、所定の試行回数毎に期待値＜Ｎｔｒｉａｌ＞と状態ｘを出力する。
期待値計算部８２は、期待値＜Ｎｔｒｉａｌ＞と状態ｘに基づいて、サンプル平均により求めたい量ｆｎ＝ｆ（ｘ（ｉ））の期待値を計算し、出力する。関数ｆ（ｘ（ｉ））は、たとえば、ｘｉ、ｘｉ×ｘｊなどである。期待値は、以下の式（１１）で表せる。
【０１６８】
【数１１】
2020135727.tif 000013
【０１６９】
式（１１）において、ｎはサンプリング回数を示し、Ｍは総サンプリング回数を示す。また、τｎは、ｎ回目のサンプリングにおける期待値＜Ｎｔｒｉａｌ＞を示す。
なお、上記の期待値計算部８２の処理を、図１に示した制御部１７が行ってもよい。
【０１７０】
図２１は、期待値計算処理を行うサンプリング装置の動作例を示すフローチャートである。なお、図１０に示した第１の実施の形態のサンプリング装置１０の動作と同じ処理については、図示及び説明を省略する。
【０１７１】
期待値計算処理が行われる場合、初期化（ステップＳ６０）の際に、期待値計算部８２は、式（１１）のｆｎτｎの積算値を示す変数Ｓと、τｎの積算値を示す変数ＴＡを０にする。また、期待値計算部８２に総サンプリング回数Ｍが設定される。さらに、サンプリング回数ｎを計数する図示しないカウンタの値が０にリセットされる。
【０１７２】
その後、サンプリングタイミングにおいて、期待値計算部８２は、サンプリング処理部８１が出力した状態ｘと期待値＜Ｎｔｒｉａｌ＞とを読み込み（ステップＳ６１）、ｆｎ＝ｆ（ｘ（ｉ））を計算する（ステップＳ６２）。その後、期待値計算部８２は、現在の変数Ｓにｆｎτｎを加算することで変数Ｓを更新し（ステップＳ６３）、現在の変数ＴＡにτｎを加算することで変数ＴＡを更新する（ステップＳ６４）。そして、期待値計算部８２は、図示しないカウンタにサンプリング回数ｎを１増加させ（ステップＳ６５）、ｎ＞Ｍであるか否かを判定する（ステップＳ６６）。ｎ＞Ｍではない場合、期待値計算部８２は、ステップＳ６１からの処理を繰り返す。
【０１７３】
ｎ＞Ｍの場合、期待値計算部８２は、Ｓ／ＴＡを計算し、式（１１）で示される期待値として出力し（ステップＳ６７）、期待値計算処理を終える。
上記のような期待値計算処理を行う期待値計算部８２は、たとえば、積和演算回路や記憶部（レジスタなど）などを用いて実現できる。また、期待値計算部８２は、たとえば、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路、ＣＰＵやＤＳＰなどのプロセッサなどであってもよい。
【０１７４】
期待値計算処理は比較的演算量が大きいため、期待値計算処理をサンプリング装置８０の外部で行う場合、サンプリング装置８０と外部装置間での状態変数ｘなどの入出力によるボトルネックが生じる可能性がある。しかし、上記のようにサンプリング装置８０内で期待値計算処理を行うことで、このボトルネックを解消できる利点がある。
【０１７５】
なお、上記各変形例は、第１乃至第３の実施の形態のサンプリング装置１０，３０，４０のそれぞれと組合せることができる。
以上、実施の形態に基づき、本発明のサンプリング装置及びサンプリング装置の制御方法の一観点について説明してきたが、これらは一例にすぎず、上記の記載に限定されるものではない。
【符号の説明】
【０１７６】
１０  サンプリング装置
１１  状態保持部
１２  エネルギー変化計算部
１３  オフセット制御部
１４  比較部
１５  フラグ計数部
１６  選択部
１７  制御部
１８  試行回数計算部
【図１】
2020135727.tif 000014
【図２】
2020135727.tif 000015
【図３】
2020135727.tif 000016
【図４】
2020135727.tif 000017
【図５】
2020135727.tif 000018
【図６】
2020135727.tif 000019
【図７】
2020135727.tif 000020
【図８】
2020135727.tif 000021
【図９】
2020135727.tif 000022
【図１０】
2020135727.tif 000023
【図１１】
2020135727.tif 000024
【図１２】
2020135727.tif 000025
【図１３】
2020135727.tif 000026
【図１４】
2020135727.tif 000027
【図１５】
2020135727.tif 000028
【図１６】
2020135727.tif 000029
【図１７】
2020135727.tif 000030
【図１８】
2020135727.tif 000031
【図１９】
2020135727.tif 000032
【図２０】
2020135727.tif 000033
【図２１】
2020135727.tif 000034