(57)【要約】【課題】マルコフ連鎖を破壊せずに効率的に評価関数の最小値を探索する。【解決手段】最適化装置１０の記憶部１１は、問題を変換した評価関数に含まれる状態変数の値を記憶し、最適化装置１０の処理部１２は、現在の状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、状態変数の値が変化することによる評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により状態変数の値を更新する処理を繰り返すことで、評価関数の最小値の探索を行う。【選択図】図１
2020205049.tif 000002
【特許請求の範囲】
【請求項１】
問題を変換した評価関数に含まれる状態変数の値を記憶する記憶部と、
現在の前記状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、前記状態変数の値が変化することによる前記評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により前記状態変数の値を更新する処理を繰り返すことで、前記評価関数の最小値の探索を行う処理部と、
を有する最適化装置。
【請求項２】
前記遷移確率分布は、前記評価関数の値の変化と逆温度との積に１を加えた値のｍ（ｍ＞１）乗の逆数で表される、請求項１に記載の最適化装置。
【請求項３】
前記遷移確率分布は、前記評価関数の値の変化と逆温度との積の２乗に１を加えた値のｍ／２（ｍ＞１）乗の逆数で表される、請求項１に記載の最適化装置。
【請求項４】
前記遷移確率分布は、前記評価関数の値の変化と逆温度との積のｍ（ｍ＞０、ただしｍ≠１）乗にマイナスを掛けた値の指数関数で表される、請求項１に記載の最適化装置。
【請求項５】
前記処理部は、前記状態変数の値を更新するたびに前記評価関数の値を計算し、計算した前記評価関数の値のうち、小さい順に複数個の値を前記記憶部に記憶させる、請求項１乃至４の何れか一項に記載の最適化装置。
【請求項６】
前記処理部は、所定のサンプリング頻度で前記評価関数の値を出力する、請求項１乃至５の何れか一項に記載の最適化装置。
【請求項７】
前記処理部は、前記状態変数の値を更新する処理の回数の増加に伴って所定のスケジュールにより０から１まで増加するアニーリング変数を用い、１から前記アニーリング変数の値を引いた値と既知の第１のハミルトニアンとの積と、前記アニーリング変数と前記評価関数の値との積との和で表される第２のハミルトニアンを、前記アニーリング変数の各値について計算して計算結果を出力する、請求項１乃至６の何れか一項に記載の最適化装置。
【請求項８】
前記処理部は、前記アニーリング変数を０から１まで増加させる複数のスケジュールの中から、前記所定のスケジュールを選択する、請求項７に記載の最適化装置。
【請求項９】
前記処理部は、複数のレプリカのそれぞれに互いに異なる温度パラメータの値を設定し、設定された前記温度パラメータの値と前記遷移確率分布に基づいた交換確率にしたがって前記複数のレプリカ間で前記温度パラメータの値を交換するレプリカ交換法により前記最小値の探索を行う、請求項１乃至８の何れか一項に記載の最適化装置。
【請求項１０】
前記処理部は、前記遷移確率分布に含まれるパラメータの値を変え、各パラメータの値における前記評価関数の値についての評価値を、前記評価関数の前記最小値の探索を行う際の計算時間よりも短い計算時間で計算し、計算結果を出力する、請求項１乃至９の何れか一項に記載の最適化装置。
【請求項１１】
最適化装置の処理部が、
記憶部に記憶されている問題を変換した評価関数に含まれる状態変数の値を取得し、
現在の前記状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、前記状態変数の値が変化することによる前記評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により前記状態変数の値を更新する処理を繰り返すことで、前記評価関数の最小値の探索を行う、
最適化装置の制御方法。
【請求項１２】
記憶部に記憶されている問題を変換した評価関数に含まれる状態変数の値を取得し、
現在の前記状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、前記状態変数の値が変化することによる前記評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により前記状態変数の値を更新する処理を繰り返すことで、前記評価関数の最小値の探索を行う、
処理をコンピュータに実行させる最適化装置の制御プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムに関する。
【背景技術】
【０００２】
自然科学や社会科学において頻出する問題として、評価関数（エネルギー関数とも呼ばれる）の最小値（または最小値を与える評価関数の状態変数の値の組合せ）を求める最小値求解問題（または組合せ最適化問題と呼ばれる）がある。近年、このような問題を磁性体のスピンの振る舞いを表すモデルであるイジングモデルで定式化する動きが加速している。この動きの技術的な基盤は、イジング型量子コンピュータの実現である。イジング型量子コンピュータは、ノイマン型コンピュータが不得意とする多変数の組合せ最適化問題を現実的な時間で解けると期待されている。一方、イジング型のコンピュータを電子回路で実装した最適化装置も開発されている（たとえば、非特許文献１参照）。
【０００３】
イジングモデルを用いた最小値求解問題の計算手法として、イジング型のエネルギー関数の最小値を、マルコフ連鎖モンテカルロ法（以下ＭＣＭＣ法、またはＭＣＭＣ計算という）を用いて探索する手法がある。ＭＣＭＣ計算では、ボルツマン分布にしたがった遷移確率でエネルギー関数の状態変数の更新（状態遷移）を行うことが一般的である。しかし、エネルギー関数には局所解となる多数の極小値が含まれ、解が一旦局所解に捕捉されると、局所解から脱出する確率が非常に低くなる。局所解からの脱出を促すため、シミュレーテッド・アニーリング法（以下ＳＡ法と略す）（たとえば、非特許文献１，２，３参照）、レプリカ交換法（拡張アンサンブル法などとも呼ばれる）（たとえば、特許文献１，２、非特許文献４参照）などの手法が知られている。
【０００４】
ＳＡ法は、遷移確率を決める式に含まれる温度パラメータの値を所定のスケジュールにしたがって、徐々に小さくしていくことで（温度を下げることに相当）、最終的にエネルギー関数の最小値を求める方法である。レプリカ交換法は、複数のレプリカのそれぞれに互いに異なる値の温度パラメータを設定し、各レプリカにおいて独立にＭＣＭＣ計算を行うとともに、所定の交換確率でレプリカ間の温度パラメータ（または状態変数）の値を交換する手法である。
【０００５】
ただ、ＳＡ法では、温度パラメータの値を非常にゆっくりと小さくしていかなければ最小値が得られないため計算の効率が悪い。また、レプリカ交換法は、各レプリカに設定される温度パラメータの値が低温領域から高温領域までくまなく動き回るようにしなければ、効率のよい最小値の探索ができないが、そのための温度パラメータの設定が難しい。一方、解が極小値に陥った場合、遷移確率を決める式に含まれるエネルギーにオフセットを加え、エネルギーが上がる方向の遷移確率を上げることで、極小値からの脱出を促進する手法（以下ダイナミックオフセット法と呼ぶ）がある。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１８−５５４１号公報
【特許文献２】特許第６４６５２３１号公報
【非特許文献】
【０００７】
【非特許文献１】Sanroku Tsukamoto, Motomu Takatsu, Satoshi Matsubara and Hirotaka Tamura, “An Accelerator Architecture for Combinatorial Optimization Problems”, FUJITSU Sci. Tech. J., Vol.53, No.5, September, 2017, pp.8-13
【非特許文献２】S. Kirkpatrick, C. D. Gelatt Jr, M. P. Vecchi, “Optimization by Simulated Annealing”, Science, Vol.220, No.4598, 13 May, 1983, pp.671-680
【非特許文献３】Constantino Tsallis, Daniel A. Stariolo, “Generalized simulated annealing”, Physica A, 233, 1996, pp.395-406
【非特許文献４】Koji Hukushima and Koji Nemoto, “Exchange Monte Carlo Method and Application to Spin Glass Simulations”, J. Phys. Soc. Jpn, Vol.65, No. 6, June, 1996, pp.1604-1608
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、極小値からの解の脱出を促進させて効率的に評価関数の最小値の探索を行うためにダイナミックオフセット法を用いる場合、エネルギーにオフセットが加えられることでマルコフ連鎖が破壊され、得られる解の精度が悪化する可能性がある。また、一旦解が極小値から脱出しても再度同じ極小値に捕捉される可能性もある。
【０００９】
１つの側面では、マルコフ連鎖を破壊せずに効率的に評価関数の最小値の探索が可能な最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
１つの実施態様では、問題を変換した評価関数に含まれる状態変数の値を記憶する記憶部と、現在の前記状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、前記状態変数の値が変化することによる前記評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により前記状態変数の値を更新する処理を繰り返すことで、前記評価関数の最小値の探索を行う処理部と、を有する最適化装置が提供される。
【００１１】
また、１つの実施態様では、最適化装置の制御方法が提供される。
また、１つの実施態様では、最適化装置の制御プログラムが提供される。
【発明の効果】
【００１２】
１つの側面では、マルコフ連鎖を破壊せずに効率的に評価関数の最小値の探索が可能となる。
【図面の簡単な説明】
【００１３】
【図１】第１の実施の形態の最適化装置の一例を示す図である。
【図２】ボルツマン分布とボルツマン分布とは異なる２つの分布の比較例を示す図である。
【図３】第２の実施の形態の最適化装置のハードウェアの一例を示す図である。
【図４】第２の実施の形態の最適化装置の機能例を示すブロック図である。
【図５】第２の実施の形態の最適化装置の制御方法の一例の処理の流れを示すフローチャートである。
【図６】情報読込処理の一例の処理の流れを示すフローチャートである。
【図７】初期化処理の一例の処理の流れを示すフローチャートである。
【図８】Ｍ−Ｈ計算処理の一例の処理の流れを示すフローチャートである。
【図９】エネルギー更新処理の一例の処理の流れを示すフローチャートである。
【図１０】アニーリング処理の一例の処理の流れを示すフローチャートである。
【図１１】ＳＡ処理の一例の処理の流れを示すフローチャートである。
【図１２】λアニール処理の一例の処理の流れを示すフローチャートである。
【図１３】レプリカ交換処理の一例の処理の流れを示すフローチャートである。
【図１４】レプリカ交換を行うペアのレプリカ番号の決定方法の一例を示す図である。
【図１５】サンプル出力処理の一例の処理の流れを示すフローチャートである。
【図１６】２種類の遷移確率分布を用いた場合のサンプリング結果の比較例を示す図である。
【図１７】温度を変更したときの２種類の遷移確率分布を用いた場合のサンプリング結果の比較例を示す図である。
【図１８】パラメータ最適化処理の一例の処理の流れを示すフローチャートである。
【図１９】評価値の出力例を示す図である。
【図２０】３種の遷移確率分布を使用した場合の定常状態の計算結果の一例を示す図である。
【図２１】探索性能の比較例を示す図である。
【図２２】ボルツマン分布を用いた場合の各レプリカにおける温度の遷移の計算結果の一例を示す図である。
【図２３】ボルツマン分布を用いた場合の各レプリカにおけるエネルギーの計算結果の一例を示す図である。
【図２４】Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合の各レプリカにおける温度の遷移の計算結果の一例を示す図である。
【図２５】Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合の各レプリカにおけるエネルギーの計算結果の一例を示す図である。
【図２６】ボルツマン分布とＰｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のレプリカ交換法の計算結果の一例を示す図である。
【図２７】図２６の計算結果において、エネルギーが小さい領域を拡大した図である。
【図２８】λアニール法の計算結果の一例を示す図である。
【図２９】デジタル回路を用いて並列試行を行う装置と、第２の実施の形態の最適化装置による計算結果の比較例を示す図である。
【図３０】連続関数であるエネルギー関数に対する探索性能の比較例を示す図である。
【発明を実施するための形態】
【００１４】
以下、発明を実施するための形態を、図面を参照しつつ説明する。
（第１の実施の形態）
図１は、第１の実施の形態の最適化装置の一例を示す図である。
【００１５】
記憶部１１は、問題を変換した評価関数（以下エネルギー関数という）に含まれる状態変数の値と状態変数の値に対応した評価関数の値（以下エネルギーという）などを記憶する。記憶部１１は、ＲＡＭ（Random Access Memory）などの揮発性の記憶装置、または、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性の記憶装置である。
【００１６】
処理部１２は、ＭＣＭＣ法により状態変数の値を更新する処理を繰り返すことでエネルギー関数の最小値の探索を行う。処理部１２は、ＣＰＵ（Central Processing Unit）、ＧＰＧＰＵ（General-Purpose computing on Graphics Processing Units）、ＤＳＰ（Digital Signal Processor）などのプロセッサである。ただし、処理部１２は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの特定用途の電子回路を含んでもよい。プロセッサは、ＲＡＭなどのメモリに記憶されたプログラムを実行する。たとえば、最適化装置１０の制御プログラムが実行される。なお、複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」ということがある。
【００１７】
最適化装置１０は、たとえば、問題を変換したイジング型のエネルギー関数の最小値（または最小値が得られる状態変数の値の組合せ）を探索するものである。
イジング型のエネルギー関数（Ｈ（｛ｘ｝））は、たとえば、以下の式（１）で定義される。
【００１８】
【数１】
2020205049.tif 000003
【００１９】
右辺の１項目は、Ｎ個の状態変数の全組合せについて、漏れと重複なく、２つの状態変数の値（０または１）と重み係数との積を積算したものである。ｘｉはｉ番目の状態変数、ｘｊはｊ番目の状態変数を表し、Ｗｉｊは、ｘｉ，ｘｊの相互作用の大きさを示す重み係数である。右辺の２項目は、各状態変数のそれぞれについてのバイアス係数（ｂｉ）と各状態変数の値との積の総和を求めたものであり、右辺の３項目（Ｃ）は定数である。
【００２０】
また、状態変数の１つであるｘｋの値が変化することによるエネルギー関数の値の変化（エネルギー差（ΔＥｋ））は、ΔＥｋ＝−（１−２ｘｋ）ｈｋと表せる。１−２ｘｋはｘｋの変化分（Δｘｋ）を表す。ｘｋが１から０に変化する場合、１−２ｘｋ＝−１であり、ｘｋが０から１に変化する場合、１−２ｘｋ＝１である。また、ｈｋはローカルフィールドと呼ばれ、以下の式（２）で表せる。
【００２１】
【数２】
2020205049.tif 000004
【００２２】
上記のようなエネルギー関数の最小値を探索する際、処理部１２は、ＭＣＭＣ計算により、複数の状態変数の何れかの値を更新する処理を繰り返す。ＭＣＭＣ計算理論の基礎となっている確率過程は以下の式（３）により記述される。
【００２３】
【数３】
2020205049.tif 000005
【００２４】
式（３）において、πｉ（｛ｘ｝，ｔ）は、｛ｘ｝（各状態変数の値を示す）で規定される状態が実現する確率である。Ｐｉ→ｊは、状態ｉから状態ｊへ遷移する遷移確率を表す。ＭＣＭＣ計算では、計算を繰り返すことで定常分布になることが保証されており、以下の式（４）によって定常状態が記述される。
【００２５】
【数４】
2020205049.tif 000006
【００２６】
ここでπｉ（｛ｘ｝，ｔ）はマルコフ連鎖が定常状態になると時間に依存しない。したがって、式（４）は一般式であるが、式（４）において、和の各項について等式が成り立つ条件がよく使用される。この条件は、詳細つり合いの原理と呼ばれ、以下の式（５）で表される。
【００２７】
【数５】
2020205049.tif 000007
【００２８】
次に、ＭＣＭＣ計算のコア部分になるメトロポリス−ヘイスティングス法（以下Ｍ−Ｈ法と略す）について説明する。
Ｍ−Ｈ法では、詳細つり合いの原理（式（５））を満たす遷移確率として、以下の式（６）が用いられる。
【００２９】
【数６】
2020205049.tif 000008
【００３０】
状態ｉから状態ｊに状態遷移するときのエネルギー差（ΔＥ＝Ｅｊ−Ｅｉ＝Ｈ（｛ｘｊ｝）−Ｈ（｛ｘｉ｝））が０より小さい場合、エネルギーが下がる状態遷移であるため、式（６）は、Ｐｉ→ｊ＝１となる。一方で、逆の状態遷移はエネルギーが上がる状態遷移であるため、Ｐｊ→ｉ＝πｉ（｛ｘ｝）／πｊ（｛ｘ｝）である。これらを式（５）に代入すると、以下の式（７）のようになり、詳細つり合いの原理を満たす。
【００３１】
【数７】
2020205049.tif 000009
【００３２】
反対にΔＥ＝Ｅｊ−Ｅｉ≧０の場合についても同様の論理で詳細つり合いの原理が満たされる。したがって、マルコフ連鎖は一意収束するので、サンプリングされる確率分布はπｉ（｛ｘ｝）になる。
【００３３】
Ｍ−Ｈ法において、π（｛ｘ｝，ｔ）の確率分布として従来のようにボルツマン分布を用いる場合、π（｛ｘ｝，ｔ）＝ｅｘｐ（−βＥ）と表せる（βは逆温度（温度パラメータの値の逆数））。そして式（６）は、以下の式（８）のように表せる。
【００３４】
【数８】
2020205049.tif 000010
【００３５】
しかしながら前述のようにボルツマン分布を用いた場合、解が一旦局所解に捕捉されると、局所解から脱出する確率が非常に低くなる。遷移確率が式（８）で表される場合、ΔＥが正に大きくなると、指数関数的にエネルギー障壁を越える確率が下がるためである。たとえば、−βΔＥ＝２０程度の場合でも、ｅｘｐ（−２０）≒２．０６×１０−９となるため、２１億回の計算回数（ＭＣＭＣ計算の回数）のうち１回程度しか、−βΔＥ＝２０となる状態遷移は生じない。
【００３６】
このため、ＳＡ法、レプリカ交換法またはダイナミックオフセット法を用いて局所解からの脱出を促進する方法があるが、前述のような技術的な問題がある。
磁性体などの物理系についての問題を扱う場合、ボルツマン分布を使わざるを得ない。なぜならば、熱平衡状態がｅｘｐ（−βＥ）の統計分布にしたがうからである。物理法則を満たさない系の動的挙動のシミュレーションに意味はないからである。また、ヘルムホルツ自由エネルギーやギブス自由エネルギーなど、系の熱力学的挙動にエントロピー効果が重要な役割を果たす場合は、エントロピーの効果を適切に取り入れるため、使用可能な遷移確率分布にも制約がかかる。
【００３７】
しかしながら、イジングモデルで表現される式（１）のようなエネルギー関数の最小値を求める問題は、単に関数の最小値求解問題とみなせばよいため、ボルツマン分布に縛られなくてもよい。
【００３８】
なお、式（１）のようなエネルギー関数において、解の個数は２Ｎ個存在する。これらをエネルギーの小さい順に並べた集合｛Ｅ｝を｛Ｅ｝＝｛Ｅ０，Ｅ１，Ｅ２，…，Ｅｉ，…，Ｅｊ，…，ＥＭ｝とする。つまり、ｉ＜ｊならば、Ｅｉ＜Ｅｊとする。ここで、ＭはＭ≦２Ｎである。なぜならば、解に重複が許されるからである。そして、Ｅ０は下に有界であるとする。つまり、負の無限大に発散せず、必ず最小値が存在するものとする。
【００３９】
以下、ボルツマン分布に代る分布を提示する前に、まず、遷移確率は以下の式（９）のように定義されるものとする。
【００４０】
【数９】
2020205049.tif 000011
【００４１】
式（９）において、Ｐｉ→ｊはΔＥ＜０のとき１になり、ΔＥ≧０のときｆ（ΔＥ）の式が適用される関数であるものとする。そしてｆ（ΔＥ）はｆ（０）＝１及び以下の式（１０）の条件を満たすものとする。
【００４２】
【数１０】
2020205049.tif 000012
【００４３】
つまり、ｆ（ΔＥ）は、現在の状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率（遷移確率）の和が１になるように規格化されている。
このとき確率行列は、以下の式（１１）のように定義される。
【００４４】
【数１１】
2020205049.tif 000013
【００４５】
また、確率行列（Ｐ）の各要素は以下の式（１２）のように定義される。
【００４６】
【数１２】
2020205049.tif 000014
【００４７】
式（１２）において、ΔＥｉｊ＝Ｅｊ−Ｅｉである。したがって、確率行列の各行について以下の式（１３）が成り立つ。
【００４８】
【数１３】
2020205049.tif 000015
【００４９】
したがって、この確率行列は固有値１をもち、定常マルコフ連鎖を形成する。
なお、上記式（９）〜式（１３）までの定式化は、イジング型の評価関数のように状態変数が０か１の値である場合（エネルギーが離散的となる場合）について行ったが、状態変数が連続変数である場合（エネルギーが連続関数となる場合）に対しても適用できる。以下にその理由を簡単に示す。連続変数が用いられる場合、エネルギーの区間［Ｅ０，Ｅ１］に対して、微小区間［Ｅｉ，Ｅｉ＋ΔＥ］における遷移確率は変わらないとして、Ｅｉにおける遷移確率が用いられる。そして、［Ｅｉ，Ｅｉ＋ΔＥ］の状態から［Ｅｊ，Ｅｊ＋ΔＥ］への遷移確率をｆ（ΔＥ）＝ｆ（Ｅｊ−Ｅｉ）とすれば、離散モデルであるイジング型の評価関数を用いた場合と、全く同じ計算となる。なお、証明の本質部分は同じであるため、状態変数が連続変数である場合についての証明を省略する。
【００５０】
なお、上記の方法は、式（５）に示した詳細つり合いの原理を一般的には破る。この方法の特殊な場合であるボルツマン分布を用いた場合（ｆ（ΔＥ）＝ｅｘｐ（−βΔＥ）とした場合）は、詳細つり合いの原理を満たす。
【００５１】
次に、ボルツマン分布に代る遷移確率分布（ｆ（ΔＥ））の例を示す。ｆ（ΔＥ）に求められる条件は、式（１０）のように規格化できるという条件であり、その条件を満たしていれば任意の遷移確率分布を適用できる。ただ、最適化装置１０の処理部１２は、最小値求解問題（最適化問題）を効率よく解くという観点で、たとえば、以下の式（１４）、式（１５）、式（１６）で表される遷移確率分布を用いる。
【００５２】
【数１４】
2020205049.tif 000016
【００５３】
【数１５】
2020205049.tif 000017
【００５４】
【数１６】
2020205049.tif 000018
【００５５】
ただし、式（１０）の規格化条件を満たすために、式（１４）〜式（１６）において、ｍ１＞１、ｍ２＞１、ｍ３＞０である。ただし、式（１４）〜式（１６）はΔＥ≧０に対して定義されるものとする。
【００５６】
なお、非特許文献３には、Ｐｑ（ｘｔ→ｘｔ＋１）＝１／［１＋（ｑ−１）βΔＥ］１／ｑ−１と表せる遷移確率分布が示されている（非特許文献３の式（５）参照）。この式は、式（１４）とは異なる。ｑ＝２のとき、式（１４）のｍ１＝１の場合と同じになるが、上記のように式（１４）では、ｍ１＞１としているためである。なお、非特許文献３では、１／［１＋（ｑ−１）βΔＥ］１／ｑ−１では、ｑ≧２の場合、規格化積分が無限大に発散する（規格化積分が存在しない）ため、確率論的に誤った遷移確率となる。また、非特許文献３では、ｑ＝１の場合、遷移確率分布はボルツマン分布となる。ｑ＜１の場合、遷移確率がべき乗で表される。たとえば、ｑ＝１／２の場合、Ｐｑ（ｘｔ→ｘｔ＋１）＝（１−βΔＥ／２）２となるが、１−βΔＥ／２＜０を満たす場合、つまりΔＥが正に大きくなると遷移確率が無限大に発散してしまう。そのため、非特許文献３ではこのような条件では、遷移確率を０にしているが、その場合、解が深い（ΔＥが大きい）極小値に捉えられた場合、脱出することが不可能になってしまう。式（１４）は、これらの状況が生じないように定式化されている。
【００５７】
なお、式（１５）で表される遷移確率分布は、ｍ２＝２の場合、コーシー分布と呼ばれる。式（１６）で表される遷移確率分布は、ｍ３＝１の場合、ボルツマン分布となり、ｍ３＝２の場合、正規分布となる。
【００５８】
したがって、式（１４）〜式（１６）のうち、式（１６）においてｍ３＝１以外の場合の分布がボルツマン分布に代る遷移確率分布となる。
なお、用いられる遷移確率分布は、エネルギー差が大きいほど、ボルツマン分布よりも遷移確率が大きくなる分布であればよいが、用途によって遷移確率分布を使い分けることが望ましい。たとえば、計算量の観点で言えば、式（１４）、及び式（１５）において、ΔＥ≫１の領域では近似的に、ｆ（ΔＥ）≒（βΔＥ）−ｍのように、ｍで定義されるべき関数になる。ｍ＝１＋δの場合、δを十分小さくとることで漸近的に以下の式（１７）のようになる。
【００５９】
【数１７】
2020205049.tif 000019
【００６０】
計算量の観点では、ボルツマン分布を遷移確率分布として用いる場合、ΔＥを超える平均の計算回数をＮＢとすると、ＮＢ＝ｅｘｐ（βΔＥ）となる。したがって、極小値から脱出するための計算回数は、ΔＥに関して指数関数的に増大する。一方、式（１４）、式（１５）を用いた場合に、ΔＥを超える平均の計算回数をＮｆとすると、Ｎｆ≒（βΔＥ）ｍであり、βΔＥに関して高々べき乗の計算回数でエネルギー関数の極小値から脱出可能である。
【００６１】
以上のように、処理部１２は、ΔＥが大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布を用いることで、ΔＥが大きい場合に遷移確率が指数関数的に下がるという問題を解決できる。
【００６２】
次に、最小値求解問題を解く際のエネルギー探索能力を効率化することを考える。
Ｍ−Ｈ法ではΔＥ＜０となる状態遷移は受諾され、エネルギーが上がる状態遷移についても確率的に受諾される。Ｍ−Ｈ法の性質上、ΔＥ≒０の領域の状態遷移が受諾されすぎるとエネルギーが上がる方向と下がる方向が均等化してしまう。その結果としてエネルギー探索能力が下がってしまう。
【００６３】
遷移確率分布としてボルツマン分布が用いられる場合、ΔＥ＝０を境にして遷移確率の微分が不連続的に変化するため、小さくエネルギーが上がる状態遷移が効率的に棄却される。これは、遷移確率分布としてボルツマン分布を使用する場合の利点である。
【００６４】
処理部１２は、ボルツマン分布のこの利点を用い、かつ、ΔＥ≫１の領域では遷移確率をΔＥのべき乗に比例するようにする。そのため、処理部１２は、式（１４）を用いる場合、ｍ１＝１＋δとして、δ＞０かつ、δ≪１とすればよい。たとえば、ｍ１＝１．０００１などとすればよい。
【００６５】
図２は、ボルツマン分布とボルツマン分布とは異なる２つの分布の比較例を示す図である。図２において横軸はｘ（＝βΔＥ）、縦軸はｆ（ｘ）（遷移確率）を表している。図２の左のグラフと右のグラフでは、ｘとｆ（ｘ）の範囲が変わっている。
【００６６】
図２では、ｆ（ｘ）＝ｅｘｐ（−ｘ）で表される分布（ボルツマン分布）と、ｆ（ｘ）＝１／（１＋ｘ２）、ｆ（ｘ）＝１／（１＋ｘ）で表される２つの分布が表されている。式（１４）に示した遷移確率分布において、ｍ１を上記のように１に近づけることで、ｆ（ｘ）＝１／（１＋ｘ）とほぼ同じ分布になる。ｆ（ｘ）＝１／（１＋ｘ２）は、式（１５）において、ｍ２＝２とした遷移確率分布に等しい。
【００６７】
図２において左側のグラフから分かるようにΔＥ≒０の領域では、ｆ（ｘ）＝１／（１＋ｘ）は、ボルツマン分布に漸近する。また、図２において右側のグラフから分かるように、ｘ（＝βΔＥ）≫１の領域においては、１／（１＋ｘ）と１／（１＋ｘ２）は、ｅｘｐ（−ｘ）よりも大きい値（遷移確率）となる。βΔＥ＝２０の場合、ボルツマン分布においては、ｆ（ｘ）は、ｅｘｐ（−２０．０）≒２．１×１０−９程度である。したがって、βΔＥ＝２０となる状態遷移は、２１億回に一度程度しか受諾されない。一方、式（１４）で表される遷移確率分布では、ｍ１を１に近づけ、１／（１＋ｘ）に近似したとき、βΔＥ＝２０の場合、遷移確率は０．０５程度であるから、βΔＥ＝２０となる状態遷移は、２０回に一度は受諾される。したがって、ボルツマン分布が用いられる場合よりも、２４００万倍も効率よく状態遷移を発生させられる。
【００６８】
図１には、最適化装置１０の制御方法の一例が示されている。
処理部１２は、問題情報を取得する（ステップＳ１）。問題情報は、たとえば、式（１）に示したエネルギー関数の重み係数（Ｗｉｊ）、バイアス係数（ｂｉ）、定数（Ｃ）などを含む。問題情報は、使用する遷移確率分布の情報（前述の式（１４）〜式（１６）のｍ１〜ｍ３の値などを含む）や、温度情報（たとえば、逆温度（β））、エネルギー関数に含まれる状態変数の値、ＭＣＭＣ計算の終了条件となる計算回数などを含んでいてもよい。問題情報は、外部から供給されてもよいし、記憶部１１に予め記憶されていてもよい。
【００６９】
次に、処理部１２は、初期化処理を行う（ステップＳ２）。初期化処理は、エネルギー関数が式（１）で表される場合、記憶部１１に記憶される状態変数であるｘ１〜ｘＮを初期化する処理を含む。ｘ１〜ｘＮは、たとえば、全て０に初期化されてもよいし、全て１に初期化されてもよい。また、ｘ１〜ｘＮは、ランダムに０と１が設定されるように初期化されてもよいし、外部から供給された値によって初期化されてもよい。また、初期化処理は、問題情報と、状態変数の初期値に基づいて、エネルギーの初期値を式（１）により計算する処理を含む。エネルギーの初期値は、現在の最小値（Ｅｍｉｎ）として記憶部１１に記憶される。
【００７０】
その後、処理部１２は更新処理を行う（ステップＳ３）。更新処理では、処理部１２は、前述の規格化条件を満たすとともにΔＥが大きいほどボルツマン分布よりも遷移確率が大きくなる、たとえば、前述の式（１４）〜式（１６）の何れかで表され遷移確率分布に基づいて、ＭＣＭＣ計算により状態変数の値を更新する。
【００７１】
たとえば、処理部１２は、ｘ１〜ｘＮの何れかをランダムに１つ選択して、選択した状態変数の値を反転（０から１または１から０に変化）させたときのΔＥを計算する。そして、処理部１２は、ΔＥが負ならば、選択した状態変数の値を反転させる（更新を確定させる）。ΔＥが正ならば、たとえば、式（１４）〜式（１６）の何れかのｆ（ΔＥ）と、０≦Ｒ≦１の一様乱数Ｒとの比較結果に基づいて、ｆ（ΔＥ）≧Ｒであるならば、処理部１２は、選択した状態変数の値を反転させる。ｆ（ΔＥ）≧Ｒでないならば、処理部１２は、選択した状態変数の値の反転を行わない。
【００７２】
また、ステップＳ３の処理において、状態変数の値の反転が生じる場合、処理部１２は、エネルギーを更新する。更新したエネルギーがこれまでの最小値である場合、処理部１２は、記憶部１１に記憶されているＥｍｉｎを更新する。なお、記憶部１１は、Ｅｍｉｎが得られたときの各状態変数の値を最適解の候補として記憶してもよい。
【００７３】
ステップＳ３の処理後、処理部１２は、ＭＣＭＣ計算の計算回数が終了条件となる所定の回数（ＮＥ）に達したか否かを判定する（ステップＳ４）。計算回数がＮＥに達していない場合、ステップＳ３からの処理が繰り返される。計算回数がＮＥに達した場合、処理部１２は、記憶部１１に記憶されている現在のＥｍｉｎを計算結果として、たとえば、外部装置（外部のコンピュータ、記憶媒体、表示装置など）に出力し（ステップＳ５）、処理を終える。
【００７４】
以上のような第１の実施の形態の最適化装置１０によれば、ΔＥが正に大きいほど、遷移確率がボルツマン分布よりも大きくなる遷移確率分布を適用するため、解が局所解から高効率で脱出可能となる。また、ダイナミックオフセット法のようにエネルギーにオフセットを加える方法ではないためマルコフ連鎖を破壊することもない。
【００７５】
以上のことから、マルコフ連鎖を破壊せずに効率的にエネルギー関数の最小値の探索が可能となる。
（第２の実施の形態）
図３は、第２の実施の形態の最適化装置のハードウェアの一例を示す図である。
【００７６】
第２の実施の形態の最適化装置２０は、たとえばコンピュータであり、ＣＰＵ２１、ＲＡＭ２２、ＨＤＤ２３、画像信号処理部２４、入力信号処理部２５、媒体リーダ２６及び通信インタフェース２７を有する。上記ユニットは、バスに接続されている。
【００７７】
ＣＰＵ２１は、プログラムの命令を実行する演算回路を含むプロセッサである。ＣＰＵ２１は、ＨＤＤ２３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ２２にロードし、プログラムを実行する。なお、ＣＰＵ２１は複数のプロセッサコアを備えてもよく、最適化装置２０は複数のプロセッサを備えてもよく、以下で説明する処理を複数のプロセッサまたはプロセッサコアを用いて並列に実行してもよい。また、複数のプロセッサの集合（マルチプロセッサ）を「プロセッサ」と呼んでもよい。
【００７８】
ＲＡＭ２２は、ＣＰＵ２１が実行するプログラムやＣＰＵ２１が演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。なお、最適化装置２０は、ＲＡＭ以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。
【００７９】
ＨＤＤ２３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、及び、データを記憶する不揮発性の記憶装置である。プログラムには、たとえば、最適化装置２０の制御プログラムが含まれる。なお、最適化装置２０は、フラッシュメモリやＳＳＤ（Solid State Drive）などの他の種類の記憶装置を備えてもよく、複数の不揮発性の記憶装置を備えてもよい。
【００８０】
画像信号処理部２４は、ＣＰＵ２１からの命令にしたがって、最適化装置２０に接続されたディスプレイ２４ａに画像を出力する。ディスプレイ２４ａとしては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、プラズマディスプレイ（ＰＤＰ：Plasma Display Panel）、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなどを用いることができる。
【００８１】
入力信号処理部２５は、最適化装置２０に接続された入力デバイス２５ａから入力信号を取得し、ＣＰＵ２１に出力する。入力デバイス２５ａとしては、マウスやタッチパネルやタッチパッドやトラックボールなどのポインティングデバイス、キーボード、リモートコントローラ、ボタンスイッチなどを用いることができる。また、最適化装置２０に、複数の種類の入力デバイスが接続されていてもよい。
【００８２】
媒体リーダ２６は、記録媒体２６ａに記録されたプログラムやデータを読み取る読み取り装置である。記録媒体２６ａとして、たとえば、磁気ディスク、光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。磁気ディスクには、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤが含まれる。光ディスクには、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）が含まれる。
【００８３】
媒体リーダ２６は、たとえば、記録媒体２６ａから読み取ったプログラムやデータを、ＲＡＭ２２やＨＤＤ２３などの他の記録媒体にコピーする。読み取られたプログラムは、たとえば、ＣＰＵ２１によって実行される。なお、記録媒体２６ａは、可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体２６ａやＨＤＤ２３を、コンピュータ読み取り可能な記録媒体ということがある。
【００８４】
通信インタフェース２７は、ネットワーク２７ａに接続され、ネットワーク２７ａを介して他の情報処理装置と通信を行うインタフェースである。通信インタフェース２７は、スイッチなどの通信装置とケーブルで接続される有線通信インタフェースでもよいし、基地局と無線リンクで接続される無線通信インタフェースでもよい。
【００８５】
次に、最適化装置２０の機能及び処理手順を説明する。
図４は、第２の実施の形態の最適化装置の機能例を示すブロック図である。
最適化装置２０は、記憶部３０、処理部３１を有する。処理部３１は、制御部３１ａ、設定読込部３１ｂ、問題設定部３１ｃ、ハミルトニアン計算部３１ｄ、エネルギー差計算部３１ｅ、乱数発生部３１ｆ、遷移確率計算部３１ｇ、ＭＣＭＣ計算実行部３１ｈ、ＳＡ実行部３１ｉ、レプリカ交換実行部３１ｊを有する。さらに、処理部３１は、エネルギー更新部３１ｋ、結果出力部３１ｌ、λアニール部３１ｍ、パラメータ最適化部３１ｎを有する。
【００８６】
なお、記憶部３０は、たとえば、ＨＤＤ２３に確保した記憶領域を用いて実装できる。処理部３１は、たとえば、ＣＰＵ２１が実行するプログラムモジュールを用いて実装できる。
【００８７】
記憶部３０は、エネルギー情報、スピン情報、レプリカ情報、温度情報、問題設定情報、ハミルトニアン情報を記憶する。
エネルギー情報は、計算されたエネルギーの初期値や、これまで計算されたエネルギーを最小値から小さい順にＮｒａｎｋ個含む。また、エネルギー情報は、Ｎｒａｎｋ個のエネルギーに対応した各状態変数の値のＮｒａｎｋ個の組合せを含んでいてもよい。スピン情報は、各状態変数の値を含む。レプリカ情報は、レプリカ交換法を実行するために用いられる情報であり、レプリカ数などを含む。温度情報は、温度パラメータの値（以下、単に温度という）または逆温度を含む。問題設定情報は、たとえば、使用する最適化方法（ＳＡ法、レプリカ交換法など）、レプリカ交換やアニーリングの頻度、サンプリング頻度、温度変更スケジュール、使用する遷移確率分布の情報、計算終了条件となるＭＣＭＣ計算の計算回数などを含む。ハミルトニアン情報は、たとえば、式（１）に示したエネルギー関数の重み係数（Ｗｉｊ）、バイアス係数（ｂｉ）、定数（Ｃ）などを含む。
【００８８】
制御部３１ａは、処理部３１の各部を制御して最小値求解処理を行う。
設定読込部３１ｂは、記憶部３０から上記の各種情報を、制御部３１ａが理解可能な形式で読み込む。
【００８９】
問題設定部３１ｃは、状態変数の初期値や、エネルギーの最小値の初期値を計算する。
ハミルトニアン計算部３１ｄは、エネルギーを計算する。
エネルギー差計算部３１ｅは、状態変数の値が反転されることによるエネルギー差を計算する。
【００９０】
乱数発生部３１ｆは、反転対象の状態変数をランダムに選択するための乱数を生成する。乱数発生部３１ｆとしてはどのようなものを使用してもよいが、ＭＣＭＣ計算の計算回数よりも十分大きい周期をもつメルセンヌツイスタなどの擬似乱数生成器などを用いることが望ましい。
【００９１】
遷移確率計算部３１ｇは、選択された遷移確率分布の式と、温度または逆温度と、計算されたエネルギー差に基づいて、遷移確率分布の値を計算する。
ＭＣＭＣ計算実行部３１ｈは、計算された遷移確率分布の値にしたがって、選択された状態変数の値を更新するか否かを判定するなどのＭＣＭＣ計算を行う。
【００９２】
ＳＡ実行部３１ｉは、ＳＡ法を実行するために、温度変更スケジュールに基づいて、温度を下げていく。
レプリカ交換実行部３１ｊは、レプリカ交換法を実行するため、後述の交換確率に基づいて、レプリカ間で温度を交換する。
【００９３】
エネルギー更新部３１ｋは、ハミルトニアン計算部３１ｄが計算したエネルギーに基づいて、記憶部３０に記憶されているエネルギー情報を更新する。
結果出力部３１ｌは、ＭＣＭＣ計算の現在の計算回数が、計算終了条件となる所定の回数（以下、ＮＥとする）に達した場合に、エネルギー情報を計算結果として出力する。なお、結果出力部３１ｌは、ＮＥよりも少ない計算回数（以下ＮＩとする）のたびに、その時点で計算されたエネルギーと、そのエネルギーに対応する各状態変数の値を出力してもよい。その場合、最適化装置２０は、サンプリング装置として機能し、ＮＩはサンプリング頻度となる。
【００９４】
λアニール部３１ｍは、ＳＡ法とは異なるアニーリング法である後述のλアニール法（パラメータアニール法、２状態遷移法とも呼ばれる）を実行する。
パラメータ最適化部３１ｎは、後述の方法により、各種パラメータを最適化する処理を行う。
【００９５】
図５は、第２の実施の形態の最適化装置の制御方法の一例の処理の流れを示すフローチャートである。
処理が開始すると、まず、設定読込部３１ｂが、記憶部３０から上記の各種情報を、制御部３１ａが理解可能な形式で読み込む（ステップＳ１０）。読み込まれる情報の例については後述する。その後、問題設定部３１ｃは、状態変数の初期化とエネルギーの初期値の計算を含む初期化処理を行う（ステップＳ１１）。そして、制御部３１ａは、たとえば、記憶部３０から読み込まれた問題設定情報にパラメータ最適化処理を行うことを指示する情報が含まれているか否かを判定する（ステップＳ１２）。制御部３１ａが、パラメータ最適化処理を行うことを指示する情報が含まれていると判定した場合、ステップＳ２４の処理が行われる。
【００９６】
制御部３１ａが、パラメータ最適化処理を行うことが指定されていないと判定した場合、制御部３１ａは、現在のＭＣＭＣ計算の計算回数（以下ステップ数と呼ぶ場合もある）を示すＮＣを１に初期化する（ステップＳ１３）。ステップＳ１３の処理後、ＭＣＭＣ計算実行部３１ｈは、ＭＣＭＣ計算により、Ｍ−Ｈ計算処理を行う（ステップＳ１４）。Ｍ−Ｈ計算処理の例については後述する。その後、エネルギー更新部３１ｋは、Ｍ−Ｈ計算処理の結果にしたがってハミルトニアン計算部３１ｄが計算したエネルギーに基づいて、記憶部３０に記憶されているエネルギー情報を更新するエネルギー更新処理を行う（ステップＳ１５）。
【００９７】
ステップＳ１５の処理後、制御部３１ａは、記憶部３０から読み込まれた問題設定情報にアニーリング（ＳＡ法またはλアニール法）を行うことを指示する情報が含まれているか否かを判定する（ステップＳ１６）。制御部３１ａが、アニーリングを行うことを指示する情報が含まれていると判定した場合、ＳＡ実行部３１ｉまたはλアニール部３１ｍによるアニーリング処理が行われる（ステップＳ１７）。アニーリング処理の例については後述する。アニーリング処理後、または、制御部３１ａが、問題設定情報にアニーリングを行うことを指示する情報が含まれていないと判定した場合、ステップＳ１８の処理が行われる。
【００９８】
ステップＳ１８の処理では、制御部３１ａは、記憶部３０から読み込まれた問題設定情報にレプリカ交換法を行うことを指示する情報が含まれているか否かを判定する（ステップＳ１８）。制御部３１ａが、レプリカ交換法を行うことを指示する情報が含まれていると判定した場合、レプリカ交換実行部３１ｊによるレプリカ交換処理が行われる（ステップＳ１９）。レプリカ交換処理の例については後述する。レプリカ交換処理後、または、制御部３１ａが、問題設定情報にレプリカ交換法を行うことを指示する情報が含まれていないと判定した場合、ステップＳ２０の処理が行われる。なお、ＳＡ法によるアニーリング処理とレプリカ交換処理とは互いに排他処理であるため、どちらか一方の選択が可能である。一方、λアニール法によるアニーリング処理とレプリカ交換処理とは同時に選択可能である。
【００９９】
ステップＳ２０の処理（サンプル出力処理）では、結果出力部３１ｌは、ＮＩ回の計算処理を行うたびに、その時点で計算されたエネルギーと、そのエネルギーに対応する各状態変数の値を出力する。サンプル出力処理の例については後述する。
【０１００】
次に、制御部３１ａは、ＮＣ＝ＮＥであるか否かを判定する（ステップＳ２１）。制御部３１ａは、ＮＣ＝ＮＥではないと判定した場合、ＮＣを＋１する（ステップＳ２２）。その後、ステップＳ１４からの処理が繰り返される。制御部３１ａがＮＣ＝ＮＥであると判定した場合、結果出力部３１ｌは、エネルギー情報を計算結果として出力する（ステップＳ２３）。その後、最適化装置２０の処理が終了する。出力されるエネルギー情報は、たとえば、これまで計算されたエネルギーを最小値から小さい順に並べたＮｒａｎｋ個のエネルギーを含む。また、エネルギー情報は、Ｎｒａｎｋ個のエネルギーに対応した各状態変数の値のＮｒａｎｋ個の組合せを含む。
【０１０１】
ステップＳ２４の処理では、パラメータ最適化部３１ｎは、後述の方法により、各種パラメータを最適化する処理を行い、その結果を表示する。その後、最適化装置２０の処理が終了する。
【０１０２】
なお、図５に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
（情報読込処理の例）
図６は、情報読込処理の一例の処理の流れを示すフローチャートである。
【０１０３】
設定読込部３１ｂは、記憶部３０から計算終了条件となる計算回数（ＮＥ）、サンプリング頻度（ＮＩ）を読み込む（ステップＳ３０，Ｓ３１）。さらに、設定読込部３１ｂは、記憶部３０から温度（たとえば、前述の式（１４）〜式（１６）のβ（＝１／Ｔ）を決めるＴ（温度（絶対温度）））を読み込む（ステップＳ３２）。また、設定読込部３１ｂは、記憶部３０から問題サイズ（状態変数の数（Ｎ））、ハミルトニアン情報（式（１）に示したエネルギー関数の重み係数（Ｗｉｊ）、バイアス係数（ｂｉ）、定数（Ｃ））を読み込む（ステップＳ３３，Ｓ３４）。
【０１０４】
また、設定読込部３１ｂは、記憶部３０からスピン初期化法（状態変数の初期値の決め方）、遷移確率分布情報を読み込む（ステップＳ３５，Ｓ３６）。遷移確率分布情報は、ＭＣＭＣ計算にどの遷移確率分布を用いるかを示す情報（分布名など）や、使用する遷移確率分布のパラメータ（たとえば、式（１４）〜式（１６）のｍ１〜ｍ３）を含む。
【０１０５】
さらに、設定読込部３１ｂは、記憶部３０から、使用する最適化方法（ＳＡ法、レプリカ交換法など）、アニーリング情報（たとえば、温度変更スケジュール、アニーリングの頻度）を読み込む（ステップＳ３７，Ｓ３８）。また、設定読込部３１ｂは、記憶部３０から、レプリカ数、レプリカ交換用の温度情報（各レプリカに設定する温度列）、アニーリングやレプリカ交換の頻度を読み込み（ステップＳ３９，Ｓ４０，Ｓ４１）、情報読込処理を終える。
【０１０６】
なお、図６に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
（初期化処理の例）
イジングモデルにおいて、スピン（状態変数）の初期化は重要である。以下に示す初期化処理の例では、問題設定部３１ｃは、４通りのスピン初期化法（“Ｅｘｔｅｒｎａｌ”、“Ｒａｎｄｏｍ”、“Ｏｎｅ”、“Ｚｅｒｏ”）の何れかを行うものとしている。
【０１０７】
図７は、初期化処理の一例の処理の流れを示すフローチャートである。
たとえば、まず問題設定部３１ｃは、スピン初期化法が“Ｅｘｔｅｒｎａｌ”であるか否かを判定する（ステップＳ５０）。問題設定部３１ｃは、スピン初期化法が“Ｅｘｔｅｒｎａｌ”であると判定した場合、最適化装置２０の外部から取得した状態変数をスピン情報の初期値として設定する（ステップＳ５１）。
【０１０８】
問題設定部３１ｃは、スピン初期化法が“Ｅｘｔｅｒｎａｌ”ではないと判定した場合、スピン初期化法が“Ｒａｎｄｏｍ”であるか否かを判定する（ステップＳ５２）。問題設定部３１ｃは、スピン初期化法が“Ｒａｎｄｏｍ”であると判定した場合、各状態変数に０または１をランダムに設定してスピン情報の初期値とする（ステップＳ５３）。たとえば、問題設定部３１ｃは、全ての状態変数を確率５０％で０に設定し、０に設定されなかった状態変数を１に設定する。
【０１０９】
問題設定部３１ｃは、スピン初期化法が“Ｒａｎｄｏｍ”ではないと判定した場合、スピン初期化法が“Ｏｎｅ”であるか否かを判定する（ステップＳ５４）。問題設定部３１ｃは、スピン初期化法が“Ｏｎｅ”であると判定した場合、全状態変数を１に設定してスピン情報の初期値とする（ステップＳ５５）。
【０１１０】
問題設定部３１ｃは、スピン初期化法が“Ｏｎｅ”ではないと判定した場合、スピン初期化法が“Ｚｅｒｏ”であると判定し、全状態変数を０に設定してスピン情報の初期値とする（ステップＳ５６）。
【０１１１】
ステップＳ５１，Ｓ５３，Ｓ５５，Ｓ５６の処理後、問題設定部３１ｃは、エネルギー（Ｅ）の初期値を計算する（ステップＳ５７）。問題設定部３１ｃは、エネルギーの初期値を、ハミルトニアン情報と、スピン情報の初期値に基づいて、式（１）により計算する。
【０１１２】
さらに、問題設定部３１ｃは、Ｅｍｉｎ＝Ｅとすることで、計算したエネルギーの初期値をエネルギーの最小値（Ｅｍｉｎ）として設定し（ステップＳ５８）、初期化処理を終える。
【０１１３】
なお、図７に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
（Ｍ−Ｈ計算処理の例）
図８は、Ｍ−Ｈ計算処理の一例の処理の流れを示すフローチャートである。
【０１１４】
乱数発生部３１ｆは、Ｎ個の状態変数の識別番号を１〜Ｎとしたとき、１〜Ｎの整数値をランダムに発生することで、反転候補の状態変数を決定する（ステップＳ６０）。
そして、エネルギー差計算部３１ｅは、反転候補の状態変数の値を反転させる場合のエネルギー差（ΔＥ）を計算する（ステップＳ６１）。反転候補の状態変数の識別番号がｋである場合、エネルギー差は、ΔＥｋ＝−（１−２ｘｋ）ｈｋである。ｈｋは、式（２）により表されるローカルフィールドである。
【０１１５】
ＭＣＭＣ計算実行部３１ｈは、エネルギー差計算部３１ｅが計算したエネルギー差が０より小さいか否かを判定する（ステップＳ６２）。ＭＣＭＣ計算実行部３１ｈがエネルギー差は０以上であると判定した場合、ＭＣＭＣ計算実行部３１ｈ（または乱数発生部３１ｆ）は、区間［０，１］で一様乱数Ｒを発生する（ステップＳ６３）。そして、ＭＣＭＣ計算実行部３１ｈは、遷移確率計算部３１ｇが計算した遷移確率分布（ｆ（ΔＥ））が、一様乱数Ｒ以上であるか否かを判定する（ステップＳ６４）。
【０１１６】
ＭＣＭＣ計算実行部３１ｈは、ｆ（ΔＥ）≧Ｒであると判定した場合、またはステップＳ６２の処理でΔＥ＜０であると判定した場合、反転候補の状態変数の値を反転させる（ステップＳ６５）。このようにして記憶部３０に記憶されるスピン情報が更新される。
【０１１７】
また、ハミルトニアン計算部３１ｄは、ステップＳ６１の処理で計算されたΔＥを用いて、エネルギーを計算（更新）する（ステップＳ６６）。
ＭＣＭＣ計算実行部３１ｈは、ステップＳ６４の処理で、ｆ（ΔＥ）＜Ｒであると判定した場合、またはステップＳ６６の処理後、１回のＭ−Ｈ計算処理を終える。
【０１１８】
なお、図８に示した処理の流れは一例である。たとえば、ステップＳ６５，Ｓ６６の処理の順序を入れ替えてもよい。
（エネルギー更新処理の一例）
図９は、エネルギー更新処理の一例の処理の流れを示すフローチャートである。
【０１１９】
エネルギー更新部３１ｋは、ハミルトニアン計算部３１ｄが計算したエネルギー（Ｅ）が、エネルギー情報として記憶部３０に記憶されているＮｒａｎｋ個のエネルギーのうち、何番目に小さいのかを示すｒａｎｋをＮｒａｎｋ＋１に初期化する（ステップＳ７０）。
【０１２０】
そして、エネルギー更新部３１ｋは、Ｎｒａｎｋ個のエネルギーの大きさの順位を示すｊを１とした後（ステップＳ７１）、Ｅ＜Ｅｍｉｎ［ｊ］であるか否かを判定する（ステップＳ７２）。Ｅｍｉｎ［ｊ］は、Ｎｒａｎｋ個のエネルギーのうち、ｊ番目に小さいエネルギーを表す。
【０１２１】
エネルギー更新部３１ｋは、Ｅ＜Ｅｍｉｎ［ｊ］であると判定した場合、ｒａｎｋ＝ｊとする（ステップＳ７３）。つまり、今回計算されたエネルギーが、Ｎｒａｎｋ個のエネルギーのうちに入る（ランキング内である）場合、その順位がｒａｎｋ＝ｊとして、たとえば、ＲＡＭ２２に保持される。
【０１２２】
エネルギー更新部３１ｋは、Ｅ≧Ｅｍｉｎ［ｊ］であると判定した場合、ｊ＝Ｎｒａｎｋであるか否かを判定し（ステップＳ７４）、ｊ＝Ｎｒａｎｋではないと判定した場合、ｊ＝ｊ＋１（ステップＳ７５）とし、ステップＳ７２からの処理を繰り返す。
【０１２３】
エネルギー更新部３１ｋは、ｒａｎｋ＝Ｎｒａｎｋ＋１のまま変わらずに、ｊ＝Ｎｒａｎｋと判定した場合（エネルギーがランキング外である場合）、エネルギー更新処理を終了する。
【０１２４】
エネルギー更新部３１ｋは、ステップＳ７３の処理後、ｒａｎｋ＝Ｎｒａｎｋであるか否かを判定する（ステップＳ７６）。エネルギー更新部３１ｋは、ｒａｎｋ＝Ｎｒａｎｋではないと判定した場合、ｊ＝ｒａｎｋ＋１とする（ステップＳ７７）。
【０１２５】
そして、エネルギー更新部３１ｋは、Ｅｍｉｎ［ｊ］＝Ｅｍｉｎ［ｊ−１］とし（ステップＳ７８）、｛ｘｊ｝＝｛ｘｊ−１｝とする（ステップＳ７９）。｛ｘｊ｝は、ｊ番目に小さいエネルギーに対応した各状態変数の値を示す。
【０１２６】
その後、エネルギー更新部３１ｋは、ｊ＝Ｎｒａｎｋであるか否かを判定し（ステップＳ８０）、ｊ＝Ｎｒａｎｋではないと判定した場合、ｊ＝ｊ＋１とし（ステップＳ８１）、ステップＳ７８からの処理を繰り返す。
【０１２７】
エネルギー更新部３１ｋは、ｊ＝Ｎｒａｎｋであると判定した場合、または、ステップＳ７６の処理で、ｒａｎｋ＝Ｎｒａｎｋであると判定した場合、Ｅｍｉｎ［ｒａｎｋ］＝Ｅとし（ステップＳ８２）、｛ｘｒａｎｋ｝＝｛ｘ｝とする（ステップＳ８３）。このようにして記憶部３０に記憶されているエネルギー情報が更新される。ステップＳ８３の処理後、エネルギー更新処理が終了する。
【０１２８】
なお、図９に示した処理の流れは一例である。適宜処理の順序が入れ替えられていてもよい。
（アニーリング処理の一例）
図１０は、アニーリング処理の一例の処理の流れを示すフローチャートである。
【０１２９】
制御部３１ａは、記憶部３０から読み込まれた問題設定情報にＳＡ法を行うことを指示する情報が含まれているか否かを判定する（ステップＳ９０）。制御部３１ａが、問題設定情報にＳＡ法を行うことを指示する情報が含まれていると判定した場合、ＳＡ実行部３１ｉは、ＳＡ処理を行う（ステップＳ９１）。制御部３１ａは、問題設定情報にＳＡ法を行うことを指示する情報が含まれていないと判定した場合、問題設定情報にλアニール法を行うことを指示する情報が含まれているか否かを判定する（ステップＳ９２）。制御部３１ａが、問題設定情報にλアニール法を行うことを指示する情報が含まれていると判定した場合、λアニール部３１ｍは、λアニール処理を行う（ステップＳ９３）。
【０１３０】
ステップＳ９１，Ｓ９３の処理後、または、ステップＳ９２の処理で問題設定情報にλアニール法を行うことを指示する情報が含まれていないと判定された場合、アニーリング処理が終了する。
【０１３１】
（ＳＡ処理の一例）
図１１は、ＳＡ処理の一例の処理の流れを示すフローチャートである。
ＳＡ実行部３１ｉは、ＭＯＤ（ＮＣ，ＮＡ）＝０であるか否かを判定する（ステップＳ１００）。ＮＣは現在の計算回数（ステップ数）であり、ＮＡは、アニーリングの頻度を示すステップ数である。ステップＳ１００の処理では、ＮＣがＮＡの倍数であるか否かが判定される。ＳＡ実行部３１ｉは、ＭＯＤ（ＮＣ，ＮＡ）＝０でないと判定した場合、ＳＡ処理を終了する。
【０１３２】
ＳＡ実行部３１ｉは、ＭＯＤ（ＮＣ，ＮＡ）＝０であると判定した場合、温度変更スケジュールを判定する（ステップＳ１０１）。
図１１の例では、４つの温度変更スケジュールの例が示されている。“ｉｎｖ”は、温度（Ｔ）を、ＮＣに反比例させて減少させる温度変更スケジュールである。“ｉｎｖｒｏｏｔ”は、Ｔの値を、ＮＣの平方根に反比例させて減少させる温度変更スケジュールである。“ｅｘｐ”は、Ｔの値を、指数型スケジュールで減少させる温度変更スケジュールである。“ｌｉｓｔ”は、Ｔの値を、リスト（任意に与えられたステップ数とＴの値との組）にしたがって変化（減少）させる温度変更スケジュールである。
【０１３３】
温度変更スケジュールが“ｉｎｖ”である場合、ＳＡ実行部３１ｉは、Ｔの値をＮＣに反比例させて減少させる（ステップＳ１０２）。温度変更スケジュールが“ｉｎｖｒｏｏｔ”である場合、ＳＡ実行部３１ｉは、Ｔの値をＮＣの平方根に反比例させて減少させる（ステップＳ１０３）。温度変更スケジュールが“ｅｘｐ”である場合、ＳＡ実行部３１ｉは、Ｔを指数型スケジュールで減少させる（ステップＳ１０４）。温度変更スケジュールが“ｌｉｓｔ”である場合、ＳＡ実行部３１ｉは、Ｔの値をリストにしたがって減少させる（ステップＳ１０２）。リストを用いることで、任意の温度変更スケジュールを設定できる。
【０１３４】
ステップＳ１０２〜Ｓ１０５の処理後、ＳＡ処理が終了する。
（λアニール処理の一例）
λアニール法では、エネルギー関数として、Ｈ（λ）＝（１−λ）ＨＡ＋λＨＢが用いられる。
【０１３５】
ＨＢは、答えを求めたいハミルトニアンであり、たとえば、式（１）が用いられる。ＨＡは初期状態のハミルトニアンであり、答えが良く知られた既知のハミルトニアンである。たとえば、１次元的なスピンの並びを考え、隣接するスピン間のみ相互作用をし、相互作用の強さは全て一定の値として与えられるものとする。こうすると、ＨＡには解析解が存在する状態を作り出せる。たとえば、ＨＡとして、全てのスピン（状態変数）の値が０のとき、エネルギー（Ｈ（λ））が最低の状態になるようなハミルトニアンを採用することができる。もちろん、全ての状態変数の値が１のとき、エネルギーが最低の状態になるようなハミルトニアンをＨＡとして採用してもよい。ＨＡの選び方は、初期状態の決め方に依存する。
【０１３６】
λ（アニーリング変数）は、０≦λ≦１であり、時間の関数としてλ＝λ（ｔ）として定義することができる。ただし、λ（ｔ）は、計算の開始時刻（ｔ＝０）にλ（０）＝０、計算の終了時刻（ｔ＝τ）にλ（τ）＝１という境界条件を満たすものとする。このため、時間の関数となるハミルトニアンは、計算の開始時刻では、Ｈ（０）＝ＨＡとなり、計算の終了時刻では、Ｈ（τ）＝ＨＢとなる。
【０１３７】
したがって、ｔ＝０にて全ての状態変数が初期化された状態から計算が開始され、計算終了時に答えの知りたいハミルトニアンの答えが得られる。これはイジング型量子コンピュータの考え方であるが、最適化装置２０は、量子論を扱わず、古典的なハミルトニアンの範囲内でエネルギーの最小値を求める手法として、λアニール法を用いる。
【０１３８】
図１２は、λアニール処理の一例の処理の流れを示すフローチャートである。
λアニール部３１ｍは、ＭＯＤ（ＮＣ，ＮＡ）＝０であるか否かを判定する（ステップＳ１１０）。ＮＣは現在の計算回数（ステップ数）であり、ＮＡは、アニーリングの頻度を示すステップ数である。ステップＳ１１０の処理では、ＮＣがＮＡの倍数であるか否かが判定されている。λアニール部３１ｍは、ＭＯＤ（ＮＣ，ＮＡ）＝０でないと判定した場合、λアニール処理を終了する。
【０１３９】
λアニール部３１ｍは、ＭＯＤ（ＮＣ，ＮＡ）＝０であると判定した場合、たとえば、読み込まれた問題設定情報に含まれるλ変更スケジュールを判定する（ステップＳ１１１）。
【０１４０】
図１２の例では、４つのλ変更スケジュールの例が示されている。“ｌｉｎｅａｒ”は、λ（ＮＣ）＝ＮＣ／ＮＥにしたがって、λ（ＮＣ）をＮＣの増加とともに増加させるλ変更スケジュールである。“ｑｕａｄｒａｔｉｃ”は、λ（ＮＣ）を（ＮＣ／ＮＥ）の２乗にしたがって、ＮＣの増加とともに増加させるλ変更スケジュールである。“ｓｑｒｔ”は、λ（ＮＣ）を（ＮＣ／ＮＥ）の平方根にしたがって、ＮＣの増加とともに増加させるλ変更スケジュールである。“ｌｉｓｔ”は、λ（ＮＣ）をｆ（ＮＣ）という関数にしたがって増加させるλ変更スケジュールである。ｆ（ＮＣ）は任意であるが、ステップ数とそのステップ数におけるλの値による組（Ｎｉ，λｉ）の列として与えられる。
【０１４１】
λ変更スケジュールが“ｌｉｎｅａｒ”である場合、λアニール部３１ｍは、λ（ＮＣ）を、λ（ＮＣ）＝ＮＣ／ＮＥを計算することで求める（ステップＳ１１２）。λ変更スケジュールが“ｑｕａｄｒａｔｉｃ”である場合、λアニール部３１ｍは、λ（ＮＣ）を、（ＮＣ／ＮＥ）の２乗を計算することで求める（ステップＳ１１３）。λ変更スケジュールが“ｓｑｒｔ”である場合、λアニール部３１ｍは、λ（ＮＣ）を、（ＮＣ／ＮＥ）の平方根を計算することで求める（ステップＳ１１４）。λ変更スケジュールが“ｌｉｓｔ”である場合、λアニール部３１ｍは、λ（ＮＣ）を、前述のｆ（ＮＣ）という関数にしたがって計算する（ステップＳ１１５）。
【０１４２】
ステップＳ１１２〜Ｓ１１５の処理後、λアニール部３１ｍは、変更前のλについてのＨ（λ）（たとえば、ハミルトニアン計算部３１ｄによって計算される）の最小値を与える各状態変数の値を記憶部３０から読み込む（ステップＳ１１６）。読み込まれた各状態変数の値は、λ変更後に再度Ｍ−Ｈ計算処理が行われるときの初期値として用いられる。ステップＳ１１６の処理後、λアニール処理が終了する。
【０１４３】
なお、λアニール処理は、以下に示すレプリカ交換処理と同時に実行させることができる。
（レプリカ交換処理の一例）
図１３は、レプリカ交換処理の一例の処理の流れを示すフローチャートである。
【０１４４】
レプリカ交換実行部３１ｊは、ＭＯＤ（ＮＣ，ＮＲ）＝０であるか否かを判定する（ステップＳ１２０）。ＮＣは現在の計算回数（ステップ数）であり、ＮＲは、レプリカ交換の頻度を示すステップ数である。ステップＳ１２０の処理では、ＮＣがＮＲの倍数であるか否かが判定されている。
【０１４５】
レプリカ交換実行部３１ｊは、ＭＯＤ（ＮＣ，ＮＲ）＝０であると判定した場合、ＭＯＤ（ＮＣ，２）＝０であるか否かを判定する（ステップＳ１２１）。ステップＳ１２１の処理では、ＮＣが偶数であるか否かが判定されている。
【０１４６】
レプリカ交換実行部３１ｊは、ＭＯＤ（ＮＣ，２）＝０であると判定した場合、温度順が偶数番目のレプリカと１つ温度が大きいレプリカによるペアを全て選択する（ステップＳ１２２）。レプリカ交換実行部３１ｊは、ＭＯＤ（ＮＣ，２）＝０でないと判定した場合、温度順が奇数番目のレプリカと１つ温度が大きいレプリカによるペアを全て選択する（ステップＳ１２３）。このようにレプリカ交換候補のレプリカペアを決めることで、レプリカ交換が効率化される。
【０１４７】
ステップＳ１２２，Ｓ１２３の処理後、レプリカ交換実行部３１ｊは、レプリカ交換を行うペアのレプリカ番号を決定する（ステップＳ１２４）。
図１４は、レプリカ交換を行うペアのレプリカ番号の決定方法の一例を示す図である。
【０１４８】
図１４において、カッコ内の上段がレプリカ番号を示し、下段がレプリカ番号のレプリカに設定されている温度を示す。ただし、Ｔ１＜Ｔ２＜Ｔ３である。レプリカ交換において交換対象のレプリカ間で温度の交換を行う場合、レプリカ番号と温度との対応関係が変化する。レプリカ交換を隣接する温度が設定されたレプリカ間で行う場合、そのようなレプリカを検出するために、ステップＳ１２４の処理では、隣接する温度（ＴｉとＴｉ＋１）に対応するレプリカ番号が検出される。
【０１４９】
図１４の、一番左側の例では、隣接する温度（Ｔ１，Ｔ２）に対応するレプリカ番号＝１，２が検出されている。レプリカ番号＝１，２のレプリカ間で温度が交換された場合、図１４の真ん中の例のようになり、隣接する温度に対応するレプリカ番号は、１，２と、１，３の２組がある。レプリカ番号＝１，２の間では同じ温度が交換されたばかりであるので、レプリカ番号＝１，３がレプリカ交換を行うペアのレプリカ番号として決定される。
【０１５０】
なお、レプリカ間で状態変数の値を交換する場合には、このような処理が不要であるが、状態変数の自由度が高くなると交換時に送受信より発生する情報量が多くなるため、上記のように温度を交換することが望ましい。
【０１５１】
その後、レプリカ交換実行部３１ｊは、レプリカ交換を行うペアのｉｄ（ｊ＝０〜Ｎｒｅｐ−１（ＮｒｅｐはステップＳ１２２，Ｓ１２３で選択されたペア数））を０にする（ステップＳ１２５）。そして、レプリカ交換実行部３１ｊは、ｉｄ＝ｊのペアに対して、レプリカ交換を実行するか否かを判定する（ステップＳ１２６）。
【０１５２】
レプリカ交換法における交換確率は、用いる遷移確率分布がボルツマン分布のように詳細つり合いの原理を満たす場合、以下の式（１８）で表せる。
【０１５３】
【数１８】
2020205049.tif 000020
【０１５４】
式（１８）において、ＰＡ（ｔ）は、レプリカ交換前の状態Ａが実現する確率である。状態Ａは、βｉが設定されているレプリカ番号＝ｉのレプリカにおける状態変数が確率分布π（βｉ，｛ｘｉ｝）にしたがい、βｊが設定されているレプリカ番号＝ｊのレプリカにおける状態変数が確率分布π（βｊ，｛ｘｊ｝）にしたがう状態である。βｉ，βｊは逆温度である。式（１８）において、ＰＢ（ｔ）は、レプリカ交換後の状態Ｂが実現する確率である。状態Ｂは、βｉが設定されているレプリカ番号＝ｉのレプリカにおける状態変数が確率分布π（βｉ，｛ｘｊ｝）にしたがい、βｊが設定されているレプリカ番号＝ｊのレプリカにおける状態変数が確率分布π（βｊ，｛ｘｉ｝）にしたがう状態である。式（１８）では、状態変数｛ｘｉ｝、｛ｘｊ｝がレプリカ番号＝ｉ，ｊのレプリカ間で交換される場合が示されているが、βｉ，βｊが交換される場合も同じ交換確率ＰＡ→Ｂとなる。
【０１５５】
一方、用いる遷移確率分布が詳細つり合いを満たさない場合、交換確率ＰＡ→Ｂは、以下の式（１９）で表される。
【０１５６】
【数１９】
2020205049.tif 000021
【０１５７】
式（１９）において、ΔＥはＥｊ−Ｅｉである。すなわち、ΔＥは、βｊが設定されているレプリカ番号＝ｊのレプリカにおけるエネルギーと、βｉが設定されているレプリカ番号＝ｉのレプリカにおけるエネルギーとの差分である。
【０１５８】
式（１９）を用いる理由を説明する。前述の式（９）により遷移確率を定義した場合、前述のように、固有値１をもつ確率行列が得られ、定常マルコフ連鎖が形成される。つまり、状態は定常状態に収束する。この定常状態に収束したときの状態密度は、遷移確率に対応した状態密度になる。いま、レプリカ交換時に全てのレプリカ同士が、式（１９）の交換確率で交換される場合、全レプリカによる系にも定常状態が存在する。レプリカ交換法などの拡張アンサンブル法を用いる場合、マルコフ連鎖が定常状態を形成することが前提であるため、式（１９）のような定常状態が得られる交換確率が用いられる。
【０１５９】
式（１９）のｆ（βｉΔＥ）として、式（１４）〜式（１６）に示したような遷移確率分布が用いられる。
ステップＳ１２６の処理では、レプリカ交換実行部３１ｊは、ｉｄ＝ｊのペアに対して、式（１９）の交換確率にしたがって、レプリカ交換を実行するか否かを判定する。たとえば、レプリカ交換実行部３１ｊは、式（１９）のｆ（βｉΔＥ）と、０≦Ｒ≦１の一様乱数Ｒとの比較結果に基づいて、ｆ（βｉΔＥ）≧Ｒであるならば、レプリカ交換実行部３１ｊは、レプリカ交換を実行すると判定する。ｆ（βｉΔＥ）≧Ｒでないならば、レプリカ交換実行部３１ｊは、レプリカ交換を実行しないと判定する。
【０１６０】
レプリカ交換実行部３１ｊは、レプリカ交換を実行すると判定した場合、ｉｄ＝ｊのペアに含まれるレプリカ間で温度を交換することで、レプリカ交換を実行する（ステップＳ１２７）。
【０１６１】
レプリカ交換実行部３１ｊは、レプリカ交換を実行しないと判定した場合、またはステップＳ１２７の処理後、ｊ＝Ｎｒｅｐ−１であるか否かを判定する（ステップＳ１２８）。レプリカ交換実行部３１ｊは、ｊ＝Ｎｒｅｐ−１であると判定した場合、レプリカ交換処理を終了し、ｊ＝Ｎｒｅｐ−１ではないと判定した場合、ｊ＝ｊ＋１とし（ステップＳ１２９）、ステップＳ１２６からの処理を繰り返す。
【０１６２】
（サンプル出力処理の一例）
図１５は、サンプル出力処理の一例の処理の流れを示すフローチャートである。図１５では、レプリカ交換処理が行われる場合の、サンプル出力処理の例が示されている。
【０１６３】
結果出力部３１ｌは、ＭＯＤ（ＮＣ，ＮＩ）＝０であるか否かを判定する（ステップＳ１３０）。ＮＣは現在の計算回数（ステップ数）であり、ＮＩは、サンプリング頻度を示すステップ数である。結果出力部３１ｌは、ＭＯＤ（ＮＣ，ＮＩ）＝０ではないと判定した場合、サンプリング出力処理を終了する。
【０１６４】
結果出力部３１ｌは、ＭＯＤ（ＮＣ，ＮＩ）＝０ではあると判定した場合、レプリカ番号（ｒ）を１に設定する（ステップＳ１３１）。ステップＳ１３１の処理後、結果出力部３１ｌは、レプリカ番号＝ｒのレプリカの現在のエネルギー（Ｅｒ）を出力し（ステップＳ１３２）、そのレプリカの現在の各状態変数の値（｛ｘｒ｝）を出力する（ステップＳ１３３）。なお、結果出力部３１ｌは、λアニール法が実行されている場合、エネルギーとして各λの値についてのＨ（λ）を出力する。
【０１６５】
さらに、結果出力部３１ｌは、ＮＣと、レプリカ番号＝ｒのレプリカに設定されている温度（Ｔｒ）と、レプリカ番号（ｒ）とを出力する（ステップＳ１３４）。
その後、結果出力部３１ｌは、ｒ＝Ｎｒであるか否かを判定する（ステップＳ１３５）。Ｎｒはレプリカ数である。結果出力部３１ｌは、ｒ＝Ｎｒであると判定した場合、サンプリング出力処理を終了し、ｒ＝Ｎｒではないと判定した場合、ｒ＝ｒ＋１とし（ステップＳ１３６）、ステップＳ１３２からの処理を繰り返す。
【０１６６】
（パラメータ最適化処理）
エネルギーの最小値をヒューリスティックに求める場合、計算者は最適な計算結果を最小計算回数または最小時間で知りたい。しかし、用いる遷移確率分布によっては解が極小値にトラップされ、簡単には脱出できない場合もある。その場合、結果としてエネルギー空間上を広くサンプリングできず（サンプリング効率が悪化し）、最小値または次善の解の求解精度と速度が悪化する。このため、最小値求解問題を計算する際には、遷移確率分布の選択が重要である。
【０１６７】
図１６は、２種類の遷移確率分布を用いた場合のサンプリング結果の比較例を示す図である。図１６の左では、式（１６）で表される遷移確率分布（以下Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布と呼ぶ場合もある）であり、ｍ３＝１／４とした場合の遷移確率分布を用いた場合のサンプリング結果が示されている。図１６の右では、式（１４）で表される遷移確率分布（以下Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布と呼ぶ場合もある）であり、ｍ１≒１とした場合の遷移確率分布を用いた場合のサンプリング結果が示されている。ただし、両分布とも、温度（絶対温度）を０．１としている。
【０１６８】
横軸はＭｏｎｔｅ  Ｃａｒｌｏ  Ｓｔｅｐｓ（前述のステップ数に相当する）を表し、縦軸はＣｏｓｔ  Ｆｕｎｃｔｉｏｎ（前述のエネルギーに相当する）を表す。
図１６の例では、サンプリング結果の分散が大きく（サンプリング範囲が広く）、よりエネルギーが低い領域を探索できるのは、Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合である。
【０１６９】
しかしながら、温度などのパラメータを変えると、結論が変わりえる。
図１７は、温度を変更したときの２種類の遷移確率分布を用いた場合のサンプリング結果の比較例を示す図である。図１７では、温度を１．０と変えた以外は、図１６と同じ条件のＴｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布と、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のサンプリング結果が示されている。
【０１７０】
温度を１．０とした場合、図１７に示されているように、よりエネルギーが低い領域を探索できるのは、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合である。一方でサンプリング結果の分散が大きいのはＴｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合である。
【０１７１】
このように、最小値求解問題を計算する際には、最適な遷移確率分布を選択することや、最適なパラメータを推定することが重要である。
一方で、最適な遷移確率分布の選択や、最適なパラメータの推定を行うための手続きは煩雑である。最小値求解問題の本計算をする前に、使用する遷移確率分布やパラメータを変えながら多数回の計算を行うことになるためである。
【０１７２】
そこで、パラメータ最適化部３１ｎは、遷移確率分布に含まれる最適なパラメータを、以下のように自動推定する機能を有する。なお、自動推定機能は、以下の性質を満たすことが望ましい。すなわち、短時間の予備計算で推定が行えること、サンプリング時のエネルギーの平均値が小さくなること、サンプリング時のエネルギーの分散が大きくなること、エネルギーの最小値がより小さくなること、という４つの性質である。
【０１７３】
図１８は、パラメータ最適化処理の一例の処理の流れを示すフローチャートである。
パラメータ最適化部３１ｎは、最適化対象の遷移確率分布の関数名を、記憶部３０から読み込む（ステップＳ１４０）。たとえば、式（１４）で表される遷移確率分布の関数名は“Ｐｏｗｅｒ  ｌａｗ”、式（１５）で表される遷移確率分布の関数名は“Ｈｙｐｅｒｂｏｌｉｃ”、式（１６）で表される遷移確率分布の関数名は“Ｔｅｍｐｅｒｅｄ−Ｇａｕｓｓｉａｎ”などと指定される。また、ユーザが指定する遷移確率分布が読み込まれる場合、“Ｕｓｅｒ−Ｄｅｆｉｎｅｄ”などと関数名が指定される。パラメータ最適化部３１ｎは、“Ｕｓｅｒ−Ｄｅｆｉｎｅｄ”を読み込む場合、（ΔＥｉ，ｆ（ΔＥｉ））の組合せも読み込み、０≦ｆ（ΔＥｉ）≦１を満たすように、各組合せによる各点を滑らかに繋いでいく。このときパラメータ最適化部３１ｎは、たとえば、キュービックスプラインなどを用いることができる。なお、パラメータ最適化部３１ｎは、ΔＥｉの最大値が与えられていて、読み込まれたΔＥｉがその最大値より大きい場合は、ｆ（ΔＥｉ）＝０とする。
【０１７４】
そして、パラメータ最適化部３１ｎは、最適化対象のパラメータ名を、記憶部３０から読み込む（ステップＳ１４１）。全ての遷移確率分布の関数において共通であるパラメータは、温度（式（１４）〜式（１６）では逆温度（β＝１／Ｔ）で表されている）であるが、式（１４）〜式（１６）の場合、ｍ１〜ｍ３も最適化対象のパラメータとして選択できる。
【０１７５】
さらに、パラメータ最適化部３１ｎは、パラメータを最適化する範囲を、記憶部３０から読み込む（ステップＳ１４２）。たとえば、パラメータが温度の場合、最小温度と最大温度が読み込まれる。なお、記憶部３０に記憶されている範囲とは異なる範囲を指定する場合には、たとえば、入力デバイス２５ａから指定すべき範囲が入力される。たとえば、パラメータが温度の場合、本実施の形態の最適化装置２０では負の絶対温度は定義されないため、範囲は０以上となる。
【０１７６】
次に、パラメータ最適化部３１ｎは、パラメータを変化させる回数（予備計算回数）Ｎｐを、記憶部３０から読み込む（ステップＳ１４３）。そして、パラメータ最適化部３１ｎは、パラメータの刻み幅ΔＰを計算する（ステップＳ１４４）。パラメータの最大値をＰｍａｘ、パラメータの最小値をＰｍｉｎとしたとき、パラメータ最適化部３１ｎは、たとえば、ΔＰ＝（Ｐｍａｘ−Ｐｍｉｎ）／（Ｎｐ−１）を計算することで、ΔＰを求める。なお、パラメータ最適化部３１ｎは、ＰｍｉｎとＰｍａｘの間を不等分に分割したΔＰを示すリストを読み込んでもよい。
【０１７７】
その後、パラメータ最適化部３１ｎは、切捨て対象の計算回数Ｎｄｕｍｐを、記憶部３０から読み込む（ステップＳ１４５）。Ｎｄｕｍｐは、数値統計の分野で“バーンイン”と呼ばれる期間を決める計算回数であり、マルコフ連鎖が定常状態に収束するまでのステップ数を意味する。計算対象の問題によってバーンインのステップ数は異なる。
【０１７８】
なお、ここまでの説明では、パラメータ最適化部３１ｎは、記憶部３０から各種情報を読み込むものとしたが、図３に示したディスプレイ２４ａに読み込み候補を表示させ、ユーザによる入力デバイス２５ａの操作によって入力された情報を読み込んでもよい。
【０１７９】
ステップＳ１４５の処理後、パラメータ最適化部３１ｎは、現在の予備計算回数を示すｉを１とし（ステップＳ１４６）、パラメータの値を計算する（ステップＳ１４７）。ステップＳ１４７の処理では、パラメータ最適化部３１ｎは、Ｐｉ＝Ｐｍｉｎ＋（ｉ−１）を計算することで、パラメータの値（Ｐｉ）を求める。
【０１８０】
その後、パラメータ最適化部３１ｎは、ＭＣＭＣ計算の現在の計算回数を示すｊを１とし（ステップＳ１４８）、ＭＣＭＣ計算によりエネルギー（Ｅｉ）とエネルギー差（ΔＥｉ）を計算する（ステップＳ１４９）。
【０１８１】
そして、パラメータ最適化部３１ｎは、ｊ＞Ｎｄｕｍｐであるか否か（計算回数がバーンイン期間を経過したか否か）を判定する（ステップＳ１５０）。パラメータ最適化部３１ｎは、ｊ＞Ｎｄｕｍｐであると判定した場合、ＥｉとΔＥｉとエネルギーの最小値を、たとえば、図３に示したＲＡＭ２３に保存する（ステップＳ１５１）。
【０１８２】
パラメータ最適化部３１ｎは、ｊ＞Ｎｄｕｍｐではないと判定した場合（マルコフ連鎖が定常状態に達していないと判定した場合）、またはステップＳ１５１の処理後、ｊ＝Ｍであるか否かを判定する（ステップＳ１５２）。Ｍは、ＭＣＭＣ計算の計算回数である。Ｍを大きくとりすぎると予備計算に時間がかかるため、バーンイン期間に相当する計算回数よりも大きく、バーンイン期間が経過するまでのデータを切り捨てた後のデータを用いて統計量を計算するためにかかる回数が指定される。モンテカルロ積分を収束するための計算量は１／（Ｍの平方根）に比例するため、推定を行える程度の計算回数が指定される。たとえば、Ｍは、１０００回から数百万回程度であり、記憶部３０から読み込まれる。これにより、以下に示す評価値の計算が、エネルギー関数の最小値の探索を行う際の計算時間（本計算時間）よりも短い時間で計算されるようにしている。
【０１８３】
パラメータ最適化部３１ｎは、ｊ＝Ｍではないと判定した場合、ｊ＝ｊ＋１とし（ステップＳ１５３）、ステップＳ１４９からの処理を繰り返す。パラメータ最適化部３１ｎは、ｊ＝Ｍであると判定した場合、保存したＥｉとΔＥｉの平均値と標準偏差を計算する（ステップＳ１５４）。なお、計算した平均値と標準偏差は、エネルギーの最小値とともに、各パラメータの値におけるエネルギー関数の値の評価値として、たとえばＲＡＭ２３に保存される。
【０１８４】
その後、パラメータ最適化部３１ｎは、ｉ＝Ｎｐであるか否かを判定する（ステップＳ１５５）。パラメータ最適化部３１ｎは、ｉ＝Ｎｐではないと判定した場合、ｉ＝ｉ＋１とし（ステップＳ１５６）、ステップＳ１４７からの処理を繰り返す。
【０１８５】
パラメータ最適化部３１ｎは、ｉ＝Ｎｐであると判定した場合、保存していた評価値を出力し（ステップＳ１５７）、パラメータ最適化処理を終える。評価値は、たとえば、図３に示したディスプレイ２４ａに出力（表示）される。なお、評価値は、記憶部３０に出力されるようにしてもよい。
【０１８６】
なお、図１８に示した処理の流れは一例である。適宜処理の順序が入れ替えられていてもよい。
図１９は、評価値の出力例を示す図である。
【０１８７】
図１９では、温度についての評価値のディスプレイ２４ａへの表示例が示されている。Ｔ＝１〜１０までＮｐ＝１０で計算を行った例が示されている。各温度についてのエネルギーの最小値（Ｅｍｉｎ）、エネルギーの平均値（Ｅａｖｅ）、エネルギーの標準偏差（Ｅｓｔｄｄｅｖ）、エネルギー差の平均値（ｄＥａｖｅ）、エネルギー差の標準偏差（ｄＥｓｔｄｄｅｖ）が評価値としてリスト化されている。また、図１９の例では、Ｅｍｉｎが小さい順に、各温度が順位付けされている。Ｎｐが大きくなると、計算結果のデータ量が増えるため、上記のように各パラメータの値に順位付けを行うことで、よいパラメータの値の判別が容易になり、解析時間を削減できる。なお、順位付けは、Ｅｍｉｎが小さい順に行うことに限定されず、Ｅａｖｅが小さい順に順位付けされるようにしてもよい。
【０１８８】
以上のような、パラメータ最適化処理を行うことで、計算者によるパラメータを決定する手間が削減され、よい計算結果が得られることが期待できるパラメータを容易に決定できる。
【０１８９】
（効果）
以上のような第２の実施の形態の最適化装置２０によれば、第１の実施の形態の最適化装置１０と同様に以下の効果が得られる。
【０１９０】
すなわち、ΔＥが（正に）大きいときの遷移確率がボルツマン分布よりも大きくなる遷移確率分布を適用するため、解が局所解から高効率で脱出可能となる。また、ダイナミックオフセット法のようにエネルギーにオフセットを加える方法ではないためマルコフ連鎖を破壊することもない。以上のことから、マルコフ連鎖を破壊せずに効率的にエネルギー関数の最小値の探索が可能となる。
【０１９１】
さらに、最適化装置２０は、レプリカ交換処理の際に、上記のような遷移確率分布を用いた交換確率を採用することで、ボルツマン分布を用いる場合よりもレプリカ交換頻度が上がる。このため、少ないレプリカ数でも温度空間を各レプリカがランダムウォークするようになる。これによりサンプリング空間が広がり、サンプリング効率が向上し、より精度の高い解を得ることができる。
【０１９２】
また、最適化装置２０は、図９に示したエネルギー更新処理において、これまで計算されたエネルギーを最小値から小さい順にＮｒａｎｋ個含むようにエネルギー情報を更新していく。これにより、エネルギー関数の真の最小値が得られない場合でも、次善の最小値（解）が得られるようになる。
【０１９３】
（計算例）
以下、上記効果を示すために、最適化装置２０による巡回セールスマン問題の計算例を示す。
【０１９４】
巡回セールスマン問題とは、Ｎ個の都市が点Ｐｉ（ｉ＝１，…，Ｎ）として定義されているとき、セールスマンが各都市を必ず１回訪問し、最後に始めの都市に戻ってくる行程の最小距離を求める問題である。ここで、点Ｐｉと点Ｐｊ間の距離ｄｉｊは、都市間の距離を表す。巡回セールスマン問題をイジングモデルに変換するために、ｉ番目に都市ａを訪問するか否かを示す状態変数をｎｉ，ａとすると、Ｎ２個の状態変数をもつイジング型のエネルギー関数（Ｅ）は以下の式（２０）のように表せる。
【０１９５】
【数２０】
2020205049.tif 000022
【０１９６】
しかし、式（２０）のままであると、セールスマンがなにも行動しない（全ての状態変数が０）という自明解や、セールスマンがｉ番目に全ての都市を訪問するような非現実解まで許容してしまうため、以下の式（２１）、式（２２）で表される拘束条件が課される。
【０１９７】
【数２１】
2020205049.tif 000023
【０１９８】
【数２２】
2020205049.tif 000024
【０１９９】
式（２１）、式（２２）で表される拘束条件をペナルティ項として、式（２０）に取り込んだ以下の式（２３）が、エネルギー最小化対象の式となる。
【０２００】
【数２３】
2020205049.tif 000025
【０２０１】
式（２３）において、ｋ１及びｋ２はペナルティの大きさを表す定数である。ペナルティの大きさをエネルギーの最小値よりも大きくしておくことで、得られる自明解や非現実解はエネルギーが大きくなり、解の候補から外れる。
【０２０２】
なお、式（２３）は、たとえば、式（１）に示した形式に変換されて計算される。ｋ１及びｋ２は重み係数（Ｗｉｊ）に反映される。
計算対象の問題として、上記のような巡回セールスマン問題を採用したのは以下の理由による。巡回セールスマン問題は、物理の問題ではないため、ボルツマン分布に縛られる必要がないという点が１つ目の理由である。また、巡回セールスマン問題は、イジングモデルに変換できるという点が２つ目の理由である。また、式（２３）において２次の項の相互作用の強さがバラバラであり、遠く離れたスピン間（状態変数間）の相互作用も強いという点が３つ目の理由である。さらに、拘束条件がペナルティ項として式（２３）に含まれるため、ペナルティ付のイジングモデル（外場ありのイジングモデル）として理想的であるという点が４つ目の理由である。
【０２０３】
図２０は、３種の遷移確率分布を使用した場合の定常状態の計算結果の一例を示す図である。横軸は、ステップ数を表し、縦軸はエネルギー（Ｅ）を表す。図２０では、遷移確率分布として、左から、ボルツマン分布、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布、Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合のサンプリング結果が示されている。
【０２０４】
使用した巡回セールスマン問題は、ＴＳＰＬＩＢのｕｌｙｓｅｅｓ１６（TSPLIB, [online], ［平成３１年４月１１日検索］、インターネット<URL: https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/>参照）である。なお、計算条件は、Ｔ＝１００．０、式（２３）のｋ１，ｋ２は１５０、ＮＥ＝１００，０００、ＮＩ＝１０、式（１４）のｍ１は１．０２、式（１６）のｍ３は０．２５である。
【０２０５】
図２０に示されているように、Ｐｏｗｅｒ  ｌａｗ型、Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合でも、ボルツマン分布を用いた場合と同様に、マルコフ連鎖を繰り返すことによって定常状態に収束していることが分かる。さらに、エネルギーの探索範囲は、ボルツマン分布を用いた場合、［１０００，４０００］程度であるが、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布では［１０００，９０００］程度まで広がる。Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合も同様である。
【０２０６】
図２１は、探索性能の比較例を示す図である。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。図２１では、遷移確率分布としてボルツマン分布を用いた場合のサンプリング結果（左上）、遷移確率分布としてボルツマン分布を用い、エネルギーにオフセットを加える手法（ダイナミックオフセット法）によるサンプリング結果（右上）が示されている。さらに、遷移確率分布として、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のサンプリング結果（左下）と、Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合のサンプリング結果（右下）が示されている。
【０２０７】
使用した巡回セールスマン問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、計算条件は、Ｔ＝０．２５、式（２３）のｋ１，ｋ２は１５０、ＮＥ＝２００，０００、ＮＩ＝１０、式（１４）のｍ１＝は１．０２、式（１６）のｍ３は０．２５である。なお、ダイナミックオフセット法の計算では、５０回連続して同じ状態に留まり続けた場合に、解が極小値に陥ったものとみなし、オフセットを加えるものとしている。
【０２０８】
図２１から明らかなように、オフセットを用いずに、ボルツマン分布を遷移確率分布とした場合、一度、解が極小値に陥った場合、解が極小値から抜け出すことはない。ダイナミックオフセット法を適用した場合、解は極小値から容易に抜け出すことができるが、解として知りたいエネルギー領域のサンプリングは数点しかない（良解が少ない）。一方、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布や、Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合、極小値から抜け出すことができるとともに、より低いエネルギー領域にも多くの解が探索できている。
【０２０９】
このため、Ｐｏｗｅｒ  ｌａｗ型やＴｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合、ボルツマン分布を用いた場合やボルツマン分布を用いたダイナミックオフセット法を適用した場合よりも、探索性能が向上することが分かった。
【０２１０】
さらに、Ｐｏｗｅｒ  ｌａｗ型や、Ｔｅｍｐｅｒｅｄ  Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合は、エネルギーにオフセットを加える手法ではないためマルコフ連鎖を壊さず、レプリカ交換法のような拡張アンサンブル法をそのまま適用できる。
【０２１１】
以下、レプリカ交換法を、遷移確率分布としてボルツマン分布を用いて実施した場合と、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いて実施した場合の計算例を示す。使用した巡回セールスマン問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、計算条件は、式（２３）のｋ１，ｋ２は１５０、ＮＥ＝２００，０００（レプリカごと）、ＮＩ＝１０、式（１４）のｍ１＝は１．０２とした。なお、レプリカ数は８つ、各レプリカに設定される温度列は、Ｔ＝１．０，３．０，５．０，１０．０，１５．０，２０．０，２５．０，３０．０であるものとした。
【０２１２】
図２２は、ボルツマン分布を用いた場合の各レプリカにおける温度の遷移の計算結果の一例を示す図である。図２２には、レプリカ番号＝１〜８までの８つのレプリカのそれぞれにおける、温度（Ｔ）の遷移の計算結果が示されている。横軸はステップ数を表し、縦軸はＴであり、正の実数を表す。
【０２１３】
図２２に示されているように、ボルツマン分布を用いた場合、各レプリカは温度空間を等確率でランダムウォークしているとは言えないことが分かる。たとえば、レプリカ番号＝５のレプリカは、ステップ数が２００付近から、Ｔ＝１．０に固定されてしまっている。
【０２１４】
図２３は、ボルツマン分布を用いた場合の各レプリカにおけるエネルギーの計算結果の一例を示す図である。図２３には、レプリカ番号＝１〜８までの８つのレプリカのそれぞれにおける、エネルギー（Ｅ）の計算結果が示されている。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。
【０２１５】
図２３に示されているように、ボルツマン分布を用いた場合、ほとんどのレプリカのエネルギーは極小値に拘束され、脱出できないことが分かる。
図２４は、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合の各レプリカにおける温度の遷移の計算結果の一例を示す図である。図２４には、レプリカ番号＝１〜８までの８つのレプリカのそれぞれにおける、温度（Ｔ）の遷移の計算結果が示されている。横軸はステップ数を表し、縦軸はＴであり、正の実数を表す。
【０２１６】
図２４に示されているように、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合、各レプリカには各温度がほぼ等確率で設定され、各レプリカは温度空間をランダムウォークしていることが分かる。
【０２１７】
図２５は、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合の各レプリカにおけるエネルギーの計算結果の一例を示す図である。図２５には、レプリカ番号＝１〜８までの８つのレプリカのそれぞれにおける、エネルギー（Ｅ）の計算結果が示されている。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。
【０２１８】
図２５に示されているように、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合、各レプリカのエネルギーは、特定の極小値に陥っても容易に脱出できていることが分かる。
図２６は、ボルツマン分布とＰｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のレプリカ交換法の計算結果の一例を示す図である。また、図２７は、図２６の計算結果において、エネルギーが小さい領域を拡大した図である。図２６、図２７において、横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。
【０２１９】
図２６から分かるように、ボルツマン分布を用いた場合（左）、エネルギーは特定の極小値に捕捉されるが、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合（右）、特定の極小値に陥っても容易に脱出できていることが分かる。また、サンプリングできるエネルギー領域も、ボルツマン分布を用いた場合よりも、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のほうが広い。
【０２２０】
さらに、図２７から分かるように、ボルツマン分布を用いた場合よりも、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のほうが得られるエネルギーの最小値も小さい。したがって、最小値探索性能も、ボルツマン分布を用いた場合よりも、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のほうが優れている。Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合、エネルギー差が大きい場合にも比較的大きな遷移確率をもつため、ボルツマン分布を用いる場合と異なり、温度の設定を粗くしても交換確率が稼げるからである。
【０２２１】
図２８は、λアニール法の計算結果の一例を示す図である。図２８の左図は、λアニール法単独の場合の計算結果を示し、右図は、λアニール法とレプリカ交換法を組合せた場合の計算結果を示す。横軸はλを表し、縦軸はエネルギー（Ｅ）（＝Ｈ（λ））を表す。
【０２２２】
使用した巡回セールスマン問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、計算条件は、式（２３）のｋ１，ｋ２＝１５０、ＮＥ＝２００，０００（レプリカ交換法の場合はレプリカごと）、ＮＩ＝１０、式（１４）のｍ１＝は１．０２とした。温度は、λアニール法単独の場合、Ｔ＝１．０、レプリカ交換法と組合せた場合、８つのレプリカに対して、Ｔ＝１．０，３．０，５．０，１０．０，１５．０，２０．０，２５．０，３０．０が設定されるものとした。
【０２２３】
λについてのアニールは、イジング型の量子コンピュータであれば、λを十分ゆっくり少しずつ増やしていけば基底状態が得られることが保証されている。しかし、量子論を扱わず、古典的なハミルトニアンの範囲内でエネルギーの最小値を求める最適化装置２０では、量子効果がないためその保証はない。図２８の左図の例では、λが大きくなるとサンプリング効率が悪化している。しかし、図２８の右図のように、λアニール法とレプリカ交換法を組合せることで、サンプリング効率を向上させることができる。これにより効率的にヒューリスティックに最小値の求解を行うことができる。
【０２２４】
図２９は、デジタル回路を用いて並列試行を行う装置と、第２の実施の形態の最適化装置による計算結果の比較例を示す図である。横軸はエネルギーの小さい順序を示すランキングを表し、縦軸はエネルギー（Ｅ）を表す。
【０２２５】
計算対象の問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、式（２３）のｋ１，ｋ２は１５０である。
計算結果４０は、たとえば、非特許文献１に示されているようなデジタル回路を用いて並列試行を行う装置を用いた場合の計算結果である（ＳＡ法とダイナミックオフセット法を適用している）。試行数（前述のＮＥに対応）は、２億回である。
【０２２６】
計算結果４１は、第２の実施の形態の最適化装置２０を用いた場合の計算結果である。なお、計算結果４１は、図５に示した処理フローに記載されているアニーリング処理やレプリカ交換処理が行われないときの計算結果を示している。ＮＥ＝２００，０００、ＮＩ＝１０である。最適化装置２０において、使用したＣＰＵ２１のクロック数は２．６０ＧＨｚであり、ＲＡＭ２２の容量は１２２８８ＭＢである。なお、プログラムはＣ＋＋で実装され、ＯｐｅｎＭＰやＭＰＩ（Message Passing Interface）による並列化は行われておらず、１ＣＰＵコアでの計算が行われた。
【０２２７】
図２９に示されているように、デジタル回路を用いて並列試行を行う装置を用いた場合よりも、最適化装置２０を用いた場合のほうが、より小さい最小値が得られ、計算精度が優れていた。また、計算時間は、前者の場合、１３．５秒、後者の場合、１．６秒であり、後者のほうが計算速度についても優れていた。
【０２２８】
次に、エネルギーが連続関数で表される場合の計算例を示す。
使用したエネルギー関数は、非特許文献３に記載されている以下の式（２４）で表される。
【０２２９】
【数２４】
2020205049.tif 000026
【０２３０】
このエネルギー関数は、Ｅ０≒５７．３２７６のときに、最小値０をとる関数である。式（２４）において、ｘｉ（状態変数）の定義域は任意の実数である。このような実数で定義されたエネルギー関数についても、式（１４）で表されるＰｏｗｅｒ  ｌａｗ型の遷移確率分布を用い、式（９）を用いたＭＣＭＣ計算を行うことで、効率的にエネルギー関数の最小値の探索が可能である。
【０２３１】
図３０は、連続関数であるエネルギー関数に対する探索性能の比較例を示す図である。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）（式（２４）のエネルギー関数の値）を表す。図３０では、遷移確率分布としてボルツマン分布を用いた場合のサンプリング結果（左）、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合のサンプリング結果（右）が示されている。
【０２３２】
計算条件は、Ｔ＝１．０、ＮＥ＝１００，０００、ＮＩ＝１０、式（１４）のｍ１＝は１．０２である。
図３０から明らかなように、ボルツマン分布を遷移確率分布とした場合、一度、解が極小値に陥った場合、解が極小値から抜け出すことはない。一方、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合、極小値から抜け出すことができることが図３０から明らかである。また、Ｐｏｗｅｒ  ｌａｗ型の遷移確率分布を用いた場合、図３０から明らかなように、一定回数同じ状態に留まるような状態を、少ない計算回数で多数求めることができる。すなわち、多数の極小値が求まるため、最小値以外にも極小値が重要な意味を持つ系では、効率的な解候補を求めることができる。そのような系の代表例として、有機化合物や生体分子などがある。
【０２３３】
（変形例）
上記の説明では、図４に示した処理部３１は、ＣＰＵ２１が実行するプログラムモジュールを用いて実装できるものとして説明したが、処理部３１の各部を、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路により実現するようにしてもよい。
【０２３４】
たとえば、ハミルトニアン計算部３１ｄは、積和演算回路を用いて式（１）を計算してもよい。ＭＣＭＣ計算実行部３１ｈについても、一様乱数と式（１４）〜式（１６）に示されているようなｆ（ΔＥ）とを比較する比較回路などを用いて実現できる。さらに、エネルギー差計算部３１ｅを、非特許文献１などのように、各状態変数の値が変化することによるエネルギー差を並列に計算する電子回路により実現するようにしてもよい。最適化装置２０は、前述のようにボルツマン分布とは異なる遷移確率分布を用いることで、高効率で解を極小値から脱出させることができるが、上記のような電子回路をアクセラレータとして用いることで、解が深い極小値に陥った場合の脱出を促進できる。同様の理由により、ＧＰＧＰＵなどをアクセラレータとして処理部３１の各部の処理を実行するために用いてもよい。
【０２３５】
なお、ＯｐｅｎＭＰやＭＰＩを用いた並列化や、ＡＶＸ（Advanced Vector eXtensions）２５６やＡＶＸ５１２など特定のＣＰＵ命令セットに対して特化したチューニングを行うことでより性能を向上できることが期待される。
【０２３６】
なお、上記では、計算対象の問題として巡回セールスマン問題を例にして説明したが、金融工学など、他の分野の問題にも適用可能である。また、ジョブスケジューリング問題もイジングモデルを用いて定式化可能であるため、各種のスケジュール最適化問題（たとえば、病院内での関係者のスケジュールを最適化する問題など）に適用することもできる。さらに、ディープラーニング分野においても、制限ボルツマンマシンはイジングモデルを用いて定式化可能であるため、制限ボルツマンマシンの最適化にも使用できる。このため、人工知能分野にも適用可能である。
【０２３７】
また、イジング型の評価関数の最小値の探索以外にも、状態変数が連続変数である評価関数の最小値の探索に対しても適用できるため、計算対象の問題も上記の分野に限定されるものではない。
【０２３８】
また、上記では評価関数の最小値を探索するものとしたが、評価関数の符号などを反転することで、評価関数の最大値を探索する手法にも拡張できる。
なお、前述のように、上記の処理内容は、最適化装置２０にプログラムを実行させることで実現できる。
【０２３９】
プログラムは、コンピュータ読み取り可能な記録媒体（たとえば、記録媒体２６ａ）に記録しておくことができる。記録媒体として、たとえば、磁気ディスク、光ディスク、光磁気ディスク、半導体メモリなどを使用できる。磁気ディスクには、ＦＤ及びＨＤＤが含まれる。光ディスクには、ＣＤ、ＣＤ−Ｒ（Recordable）／ＲＷ（Rewritable）、ＤＶＤ及びＤＶＤ−Ｒ／ＲＷが含まれる。プログラムは、可搬型の記録媒体に記録されて配布されることがある。その場合、可搬型の記録媒体から他の記録媒体（たとえば、ＨＤＤ２３）にプログラムをコピーして実行してもよい。
【０２４０】
以上、実施の形態に基づき、本発明の最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムの一観点について説明してきたが、これらは一例にすぎず、上記の記載に限定されるものではない。
【符号の説明】
【０２４１】
１０  最適化装置
１１  記憶部
１２  処理部
【図１】
2020205049.tif 000027
【図２】
2020205049.tif 000028
【図３】
2020205049.tif 000029
【図４】
2020205049.tif 000030
【図５】
2020205049.tif 000031
【図６】
2020205049.tif 000032
【図７】
2020205049.tif 000033
【図８】
2020205049.tif 000034
【図９】
2020205049.tif 000035
【図１０】
2020205049.tif 000036
【図１１】
2020205049.tif 000037
【図１２】
2020205049.tif 000038
【図１３】
2020205049.tif 000039
【図１４】
2020205049.tif 000040
【図１５】
2020205049.tif 000041
【図１６】
2020205049.tif 000042
【図１７】
2020205049.tif 000043
【図１８】
2020205049.tif 000044
【図１９】
2020205049.tif 000045
【図２０】
2020205049.tif 000046
【図２１】
2020205049.tif 000047
【図２２】
2020205049.tif 000048
【図２３】
2020205049.tif 000049
【図２４】
2020205049.tif 000050
【図２５】
2020205049.tif 000051
【図２６】
2020205049.tif 000052
【図２７】
2020205049.tif 000053
【図２８】
2020205049.tif 000054
【図２９】
2020205049.tif 000055
【図３０】
2020205049.tif 000056