JP7488457B2 - 最適化装置、最適化装置の制御方法及び最適化装置の制御プログラム - Google Patents

Description

本発明は、最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムに関する。

自然科学や社会科学において頻出する問題として、評価関数（エネルギー関数とも呼ばれる）の最小値（または最小値を与える評価関数の状態変数の値の組合せ）を求める最小値求解問題（または組合せ最適化問題と呼ばれる）がある。近年、このような問題を磁性体のスピンの振る舞いを表すモデルであるイジングモデルで定式化する動きが加速している。この動きの技術的な基盤は、イジング型量子コンピュータの実現である。イジング型量子コンピュータは、ノイマン型コンピュータが不得意とする多変数の組合せ最適化問題を現実的な時間で解けると期待されている。一方、イジング型のコンピュータを電子回路で実装した最適化装置も開発されている（たとえば、非特許文献１参照）。

イジングモデルを用いた最小値求解問題の計算手法として、イジング型のエネルギー関数の最小値を、マルコフ連鎖モンテカルロ法（以下ＭＣＭＣ法、またはＭＣＭＣ計算という）を用いて探索する手法がある。ＭＣＭＣ計算では、ボルツマン分布にしたがった遷移確率でエネルギー関数の状態変数の更新（状態遷移）を行うことが一般的である。しかし、エネルギー関数には局所解となる多数の極小値が含まれ、解が一旦局所解に捕捉されると、局所解から脱出する確率が非常に低くなる。局所解からの脱出を促すため、シミュレーテッド・アニーリング法（以下ＳＡ法と略す）（たとえば、非特許文献１，２，３参照）、レプリカ交換法（拡張アンサンブル法などとも呼ばれる）（たとえば、特許文献１，２、非特許文献４参照）などの手法が知られている。

ＳＡ法は、遷移確率を決める式に含まれる温度パラメータの値を所定のスケジュールにしたがって、徐々に小さくしていくことで（温度を下げることに相当）、最終的にエネルギー関数の最小値を求める方法である。レプリカ交換法は、複数のレプリカのそれぞれに互いに異なる値の温度パラメータを設定し、各レプリカにおいて独立にＭＣＭＣ計算を行うとともに、所定の交換確率でレプリカ間の温度パラメータ（または状態変数）の値を交換する手法である。

ただ、ＳＡ法では、温度パラメータの値を非常にゆっくりと小さくしていかなければ最小値が得られないため計算の効率が悪い。また、レプリカ交換法は、各レプリカに設定される温度パラメータの値が低温領域から高温領域までくまなく動き回るようにしなければ、効率のよい最小値の探索ができないが、そのための温度パラメータの設定が難しい。一方、解が極小値に陥った場合、遷移確率を決める式に含まれるエネルギーにオフセットを加え、エネルギーが上がる方向の遷移確率を上げることで、極小値からの脱出を促進する手法（以下ダイナミックオフセット法と呼ぶ）がある。

特開２０１８－５５４１号公報 特許第６４６５２３１号公報

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）などの特定用途の電子回路を含んでもよい。プロセッサは、ＲＡＭなどのメモリに記憶されたプログラムを実行する。たとえば、最適化装置１０の制御プログラムが実行される。なお、複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」ということがある。

最適化装置１０は、たとえば、問題を変換したイジング型のエネルギー関数の最小値（または最小値が得られる状態変数の値の組合せ）を探索するものである。
イジング型のエネルギー関数（Ｈ（｛ｘ｝））は、たとえば、以下の式（１）で定義される。

右辺の１項目は、Ｎ個の状態変数の全組合せについて、漏れと重複なく、２つの状態変数の値（０または１）と重み係数との積を積算したものである。ｘ_ｉはｉ番目の状態変数、ｘ_ｊはｊ番目の状態変数を表し、Ｗ_ｉｊは、ｘ_ｉ，ｘ_ｊの相互作用の大きさを示す重み係数である。右辺の２項目は、各状態変数のそれぞれについてのバイアス係数（ｂ_ｉ）と各状態変数の値との積の総和を求めたものであり、右辺の３項目（Ｃ）は定数である。

また、状態変数の１つであるｘ_ｋの値が変化することによるエネルギー関数の値の変化（エネルギー差（ΔＥ_ｋ））は、ΔＥ_ｋ＝－（１－２ｘ_ｋ）ｈ_ｋと表せる。１－２ｘ_ｋはｘ_ｋの変化分（Δｘ_ｋ）を表す。ｘ_ｋが１から０に変化する場合、１－２ｘ_ｋ＝－１であり、ｘ_ｋが０から１に変化する場合、１－２ｘ_ｋ＝１である。また、ｈ_ｋはローカルフィールドと呼ばれ、以下の式（２）で表せる。

上記のようなエネルギー関数の最小値を探索する際、処理部１２は、ＭＣＭＣ計算により、複数の状態変数の何れかの値を更新する処理を繰り返す。ＭＣＭＣ計算理論の基礎となっている確率過程は以下の式（３）により記述される。

式（３）において、π_ｉ（｛ｘ｝，ｔ）は、｛ｘ｝（各状態変数の値を示す）で規定される状態が実現する確率である。Ｐ_ｉ→ｊは、状態ｉから状態ｊへ遷移する遷移確率を表す。ＭＣＭＣ計算では、計算を繰り返すことで定常分布になることが保証されており、以下の式（４）によって定常状態が記述される。

ここでπ_ｉ（｛ｘ｝，ｔ）はマルコフ連鎖が定常状態になると時間に依存しない。したがって、式（４）は一般式であるが、式（４）において、和の各項について等式が成り立つ条件がよく使用される。この条件は、詳細つり合いの原理と呼ばれ、以下の式（５）で表される。

次に、ＭＣＭＣ計算のコア部分になるメトロポリス－ヘイスティングス法（以下Ｍ－Ｈ法と略す）について説明する。
Ｍ－Ｈ法では、詳細つり合いの原理（式（５））を満たす遷移確率として、以下の式（６）が用いられる。

状態ｉから状態ｊに状態遷移するときのエネルギー差（ΔＥ＝Ｅ_ｊ－Ｅ_ｉ＝Ｈ（｛ｘ_ｊ｝）－Ｈ（｛ｘ_ｉ｝））が０より小さい場合、エネルギーが下がる状態遷移であるため、式（６）は、Ｐ_ｉ→ｊ＝１となる。一方で、逆の状態遷移はエネルギーが上がる状態遷移であるため、Ｐ_ｊ→ｉ＝π_ｉ（｛ｘ｝）／π_ｊ（｛ｘ｝）である。これらを式（５）に代入すると、以下の式（７）のようになり、詳細つり合いの原理を満たす。

反対にΔＥ＝Ｅ_ｊ－Ｅ_ｉ≧０の場合についても同様の論理で詳細つり合いの原理が満たされる。したがって、マルコフ連鎖は一意収束するので、サンプリングされる確率分布はπ_ｉ（｛ｘ｝）になる。

Ｍ－Ｈ法において、π（｛ｘ｝，ｔ）の確率分布として従来のようにボルツマン分布を用いる場合、π（｛ｘ｝，ｔ）＝ｅｘｐ（－βＥ）と表せる（βは逆温度（温度パラメータの値の逆数））。そして式（６）は、以下の式（８）のように表せる。

しかしながら前述のようにボルツマン分布を用いた場合、解が一旦局所解に捕捉されると、局所解から脱出する確率が非常に低くなる。遷移確率が式（８）で表される場合、ΔＥが正に大きくなると、指数関数的にエネルギー障壁を越える確率が下がるためである。たとえば、－βΔＥ＝２０程度の場合でも、ｅｘｐ（－２０）≒２．０６×１０^－９となるため、２１億回の計算回数（ＭＣＭＣ計算の回数）のうち１回程度しか、－βΔＥ＝２０となる状態遷移は生じない。

このため、ＳＡ法、レプリカ交換法またはダイナミックオフセット法を用いて局所解からの脱出を促進する方法があるが、前述のような技術的な問題がある。
磁性体などの物理系についての問題を扱う場合、ボルツマン分布を使わざるを得ない。なぜならば、熱平衡状態がｅｘｐ（－βＥ）の統計分布にしたがうからである。物理法則を満たさない系の動的挙動のシミュレーションに意味はないからである。また、ヘルムホルツ自由エネルギーやギブス自由エネルギーなど、系の熱力学的挙動にエントロピー効果が重要な役割を果たす場合は、エントロピーの効果を適切に取り入れるため、使用可能な遷移確率分布にも制約がかかる。

しかしながら、イジングモデルで表現される式（１）のようなエネルギー関数の最小値を求める問題は、単に関数の最小値求解問題とみなせばよいため、ボルツマン分布に縛られなくてもよい。

なお、式（１）のようなエネルギー関数において、解の個数は２^Ｎ個存在する。これらをエネルギーの小さい順に並べた集合｛Ｅ｝を｛Ｅ｝＝｛Ｅ_０，Ｅ_１，Ｅ_２，…，Ｅ_ｉ，…，Ｅ_ｊ，…，Ｅ_Ｍ｝とする。つまり、ｉ＜ｊならば、Ｅ_ｉ＜Ｅ_ｊとする。ここで、ＭはＭ≦２^Ｎである。なぜならば、解に重複が許されるからである。そして、Ｅ_０は下に有界であるとする。つまり、負の無限大に発散せず、必ず最小値が存在するものとする。

以下、ボルツマン分布に代る分布を提示する前に、まず、遷移確率は以下の式（９）のように定義されるものとする。

式（９）において、Ｐ_ｉ→ｊはΔＥ＜０のとき１になり、ΔＥ≧０のときｆ（ΔＥ）の式が適用される関数であるものとする。そしてｆ（ΔＥ）はｆ（０）＝１及び以下の式（１０）の条件を満たすものとする。

つまり、ｆ（ΔＥ）は、現在の状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率（遷移確率）の和が１になるように規格化されている。
このとき確率行列は、以下の式（１１）のように定義される。

また、確率行列（Ｐ）の各要素は以下の式（１２）のように定義される。

式（１２）において、ΔＥ_ｉｊ＝Ｅ_ｊ－Ｅ_ｉである。したがって、確率行列の各行について以下の式（１３）が成り立つ。

したがって、この確率行列は固有値１をもち、定常マルコフ連鎖を形成する。
なお、上記式（９）～式（１３）までの定式化は、イジング型の評価関数のように状態変数が０か１の値である場合（エネルギーが離散的となる場合）について行ったが、状態変数が連続変数である場合（エネルギーが連続関数となる場合）に対しても適用できる。以下にその理由を簡単に示す。連続変数が用いられる場合、エネルギーの区間［Ｅ_０，Ｅ_１］に対して、微小区間［Ｅ_ｉ，Ｅ_ｉ＋ΔＥ］における遷移確率は変わらないとして、Ｅ_ｉにおける遷移確率が用いられる。そして、［Ｅ_ｉ，Ｅ_ｉ＋ΔＥ］の状態から［Ｅ_ｊ，Ｅ_ｊ＋ΔＥ］への遷移確率をｆ（ΔＥ）＝ｆ（Ｅ_ｊ－Ｅ_ｉ）とすれば、離散モデルであるイジング型の評価関数を用いた場合と、全く同じ計算となる。なお、証明の本質部分は同じであるため、状態変数が連続変数である場合についての証明を省略する。

なお、上記の方法は、式（５）に示した詳細つり合いの原理を一般的には破る。この方法の特殊な場合であるボルツマン分布を用いた場合（ｆ（ΔＥ）＝ｅｘｐ（－βΔＥ）とした場合）は、詳細つり合いの原理を満たす。

次に、ボルツマン分布に代る遷移確率分布（ｆ（ΔＥ））の例を示す。ｆ（ΔＥ）に求められる条件は、式（１０）のように規格化できるという条件であり、その条件を満たしていれば任意の遷移確率分布を適用できる。ただ、最適化装置１０の処理部１２は、最小値求解問題（最適化問題）を効率よく解くという観点で、たとえば、以下の式（１４）、式（１５）、式（１６）で表される遷移確率分布を用いる。

ただし、式（１０）の規格化条件を満たすために、式（１４）～式（１６）において、ｍ_１＞１、ｍ_２＞１、ｍ_３＞０である。ただし、式（１４）～式（１６）はΔＥ≧０に対して定義されるものとする。

なお、非特許文献３には、Ｐ_ｑ（ｘ_ｔ→ｘ_ｔ＋１）＝１／［１＋（ｑ－１）βΔＥ］^{１／ｑ－１}と表せる遷移確率分布が示されている（非特許文献３の式（５）参照）。この式は、式（１４）とは異なる。ｑ＝２のとき、式（１４）のｍ_１＝１の場合と同じになるが、上記のように式（１４）では、ｍ_１＞１としているためである。なお、非特許文献３では、１／［１＋（ｑ－１）βΔＥ］^{１／ｑ－１}では、ｑ≧２の場合、規格化積分が無限大に発散する（規格化積分が存在しない）ため、確率論的に誤った遷移確率となる。また、非特許文献３では、ｑ＝１の場合、遷移確率分布はボルツマン分布となる。ｑ＜１の場合、遷移確率がべき乗で表される。たとえば、ｑ＝１／２の場合、Ｐ_ｑ（ｘ_ｔ→ｘ_ｔ＋１）＝（１－βΔＥ／２）^２となるが、１－βΔＥ／２＜０を満たす場合、つまりΔＥが正に大きくなると遷移確率が無限大に発散してしまう。そのため、非特許文献３ではこのような条件では、遷移確率を０にしているが、その場合、解が深い（ΔＥが大きい）極小値に捉えられた場合、脱出することが不可能になってしまう。式（１４）は、これらの状況が生じないように定式化されている。

なお、式（１５）で表される遷移確率分布は、ｍ_２＝２の場合、コーシー分布と呼ばれる。式（１６）で表される遷移確率分布は、ｍ_３＝１の場合、ボルツマン分布となり、ｍ_３＝２の場合、正規分布となる。

したがって、式（１４）～式（１６）のうち、式（１６）においてｍ_３＝１以外の場合の分布がボルツマン分布に代る遷移確率分布となる。
なお、用いられる遷移確率分布は、エネルギー差が大きいほど、ボルツマン分布よりも遷移確率が大きくなる分布であればよいが、用途によって遷移確率分布を使い分けることが望ましい。たとえば、計算量の観点で言えば、式（１４）、及び式（１５）において、ΔＥ≫１の領域では近似的に、ｆ（ΔＥ）≒（βΔＥ）^－ｍのように、ｍで定義されるべき関数になる。ｍ＝１＋δの場合、δを十分小さくとることで漸近的に以下の式（１７）のようになる。

計算量の観点では、ボルツマン分布を遷移確率分布として用いる場合、ΔＥを超える平均の計算回数をＮ_Ｂとすると、Ｎ_Ｂ＝ｅｘｐ（βΔＥ）となる。したがって、極小値から脱出するための計算回数は、ΔＥに関して指数関数的に増大する。一方、式（１４）、式（１５）を用いた場合に、ΔＥを超える平均の計算回数をＮ_ｆとすると、Ｎ_ｆ≒（βΔＥ）^ｍであり、βΔＥに関して高々べき乗の計算回数でエネルギー関数の極小値から脱出可能である。

以上のように、処理部１２は、ΔＥが大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布を用いることで、ΔＥが大きい場合に遷移確率が指数関数的に下がるという問題を解決できる。

次に、最小値求解問題を解く際のエネルギー探索能力を効率化することを考える。
Ｍ－Ｈ法ではΔＥ＜０となる状態遷移は受諾され、エネルギーが上がる状態遷移についても確率的に受諾される。Ｍ－Ｈ法の性質上、ΔＥ≒０の領域の状態遷移が受諾されすぎるとエネルギーが上がる方向と下がる方向が均等化してしまう。その結果としてエネルギー探索能力が下がってしまう。

遷移確率分布としてボルツマン分布が用いられる場合、ΔＥ＝０を境にして遷移確率の微分が不連続的に変化するため、小さくエネルギーが上がる状態遷移が効率的に棄却される。これは、遷移確率分布としてボルツマン分布を使用する場合の利点である。

処理部１２は、ボルツマン分布のこの利点を用い、かつ、ΔＥ≫１の領域では遷移確率をΔＥのべき乗に比例するようにする。そのため、処理部１２は、式（１４）を用いる場合、ｍ_１＝１＋δとして、δ＞０かつ、δ≪１とすればよい。たとえば、ｍ_１＝１．０００１などとすればよい。

図２は、ボルツマン分布とボルツマン分布とは異なる２つの分布の比較例を示す図である。図２において横軸はｘ（＝βΔＥ）、縦軸はｆ（ｘ）（遷移確率）を表している。図２の左のグラフと右のグラフでは、ｘとｆ（ｘ）の範囲が変わっている。

図２では、ｆ（ｘ）＝ｅｘｐ（－ｘ）で表される分布（ボルツマン分布）と、ｆ（ｘ）＝１／（１＋ｘ^２）、ｆ（ｘ）＝１／（１＋ｘ）で表される２つの分布が表されている。式（１４）に示した遷移確率分布において、ｍ_１を上記のように１に近づけることで、ｆ（ｘ）＝１／（１＋ｘ）とほぼ同じ分布になる。ｆ（ｘ）＝１／（１＋ｘ^２）は、式（１５）において、ｍ_２＝２とした遷移確率分布に等しい。

図２において左側のグラフから分かるようにΔＥ≒０の領域では、ｆ（ｘ）＝１／（１＋ｘ）は、ボルツマン分布に漸近する。また、図２において右側のグラフから分かるように、ｘ（＝βΔＥ）≫１の領域においては、１／（１＋ｘ）と１／（１＋ｘ^２）は、ｅｘｐ（－ｘ）よりも大きい値（遷移確率）となる。βΔＥ＝２０の場合、ボルツマン分布においては、ｆ（ｘ）は、ｅｘｐ（－２０．０）≒２．１×１０^－９程度である。したがって、βΔＥ＝２０となる状態遷移は、２１億回に一度程度しか受諾されない。一方、式（１４）で表される遷移確率分布では、ｍ_１を１に近づけ、１／（１＋ｘ）に近似したとき、βΔＥ＝２０の場合、遷移確率は０．０５程度であるから、βΔＥ＝２０となる状態遷移は、２０回に一度は受諾される。したがって、ボルツマン分布が用いられる場合よりも、２４００万倍も効率よく状態遷移を発生させられる。

図１には、最適化装置１０の制御方法の一例が示されている。
処理部１２は、問題情報を取得する（ステップＳ１）。問題情報は、たとえば、式（１）に示したエネルギー関数の重み係数（Ｗ_ｉｊ）、バイアス係数（ｂ_ｉ）、定数（Ｃ）などを含む。問題情報は、使用する遷移確率分布の情報（前述の式（１４）～式（１６）のｍ_１～ｍ_３の値などを含む）や、温度情報（たとえば、逆温度（β））、エネルギー関数に含まれる状態変数の値、ＭＣＭＣ計算の終了条件となる計算回数などを含んでいてもよい。問題情報は、外部から供給されてもよいし、記憶部１１に予め記憶されていてもよい。

次に、処理部１２は、初期化処理を行う（ステップＳ２）。初期化処理は、エネルギー関数が式（１）で表される場合、記憶部１１に記憶される状態変数であるｘ_１～ｘ_Ｎを初期化する処理を含む。ｘ_１～ｘ_Ｎは、たとえば、全て０に初期化されてもよいし、全て１に初期化されてもよい。また、ｘ_１～ｘ_Ｎは、ランダムに０と１が設定されるように初期化されてもよいし、外部から供給された値によって初期化されてもよい。また、初期化処理は、問題情報と、状態変数の初期値に基づいて、エネルギーの初期値を式（１）により計算する処理を含む。エネルギーの初期値は、現在の最小値（Ｅ_ｍｉｎ）として記憶部１１に記憶される。

その後、処理部１２は更新処理を行う（ステップＳ３）。更新処理では、処理部１２は、前述の規格化条件を満たすとともにΔＥが大きいほどボルツマン分布よりも遷移確率が大きくなる、たとえば、前述の式（１４）～式（１６）の何れかで表され遷移確率分布に基づいて、ＭＣＭＣ計算により状態変数の値を更新する。

たとえば、処理部１２は、ｘ_１～ｘ_Ｎの何れかをランダムに１つ選択して、選択した状態変数の値を反転（０から１または１から０に変化）させたときのΔＥを計算する。そして、処理部１２は、ΔＥが負ならば、選択した状態変数の値を反転させる（更新を確定させる）。ΔＥが正ならば、たとえば、式（１４）～式（１６）の何れかのｆ（ΔＥ）と、０≦Ｒ≦１の一様乱数Ｒとの比較結果に基づいて、ｆ（ΔＥ）≧Ｒであるならば、処理部１２は、選択した状態変数の値を反転させる。ｆ（ΔＥ）≧Ｒでないならば、処理部１２は、選択した状態変数の値の反転を行わない。

また、ステップＳ３の処理において、状態変数の値の反転が生じる場合、処理部１２は、エネルギーを更新する。更新したエネルギーがこれまでの最小値である場合、処理部１２は、記憶部１１に記憶されているＥ_ｍｉｎを更新する。なお、記憶部１１は、Ｅ_ｍｉｎが得られたときの各状態変数の値を最適解の候補として記憶してもよい。

ステップＳ３の処理後、処理部１２は、ＭＣＭＣ計算の計算回数が終了条件となる所定の回数（Ｎ_Ｅ）に達したか否かを判定する（ステップＳ４）。計算回数がＮ_Ｅに達していない場合、ステップＳ３からの処理が繰り返される。計算回数がＮ_Ｅに達した場合、処理部１２は、記憶部１１に記憶されている現在のＥ_ｍｉｎを計算結果として、たとえば、外部装置（外部のコンピュータ、記憶媒体、表示装置など）に出力し（ステップＳ５）、処理を終える。

以上のような第１の実施の形態の最適化装置１０によれば、ΔＥが正に大きいほど、遷移確率がボルツマン分布よりも大きくなる遷移確率分布を適用するため、解が局所解から高効率で脱出可能となる。また、ダイナミックオフセット法のようにエネルギーにオフセットを加える方法ではないためマルコフ連鎖を破壊することもない。

以上のことから、マルコフ連鎖を破壊せずに効率的にエネルギー関数の最小値の探索が可能となる。
（第２の実施の形態）
図３は、第２の実施の形態の最適化装置のハードウェアの一例を示す図である。

第２の実施の形態の最適化装置２０は、たとえばコンピュータであり、ＣＰＵ２１、ＲＡＭ２２、ＨＤＤ２３、画像信号処理部２４、入力信号処理部２５、媒体リーダ２６及び通信インタフェース２７を有する。上記ユニットは、バスに接続されている。

ＣＰＵ２１は、プログラムの命令を実行する演算回路を含むプロセッサである。ＣＰＵ２１は、ＨＤＤ２３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ２２にロードし、プログラムを実行する。なお、ＣＰＵ２１は複数のプロセッサコアを備えてもよく、最適化装置２０は複数のプロセッサを備えてもよく、以下で説明する処理を複数のプロセッサまたはプロセッサコアを用いて並列に実行してもよい。また、複数のプロセッサの集合（マルチプロセッサ）を「プロセッサ」と呼んでもよい。

ＲＡＭ２２は、ＣＰＵ２１が実行するプログラムやＣＰＵ２１が演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。なお、最適化装置２０は、ＲＡＭ以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。

ＨＤＤ２３は、ＯＳ（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）が含まれる。

媒体リーダ２６は、たとえば、記録媒体２６ａから読み取ったプログラムやデータを、ＲＡＭ２２やＨＤＤ２３などの他の記録媒体にコピーする。読み取られたプログラムは、たとえば、ＣＰＵ２１によって実行される。なお、記録媒体２６ａは、可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体２６ａやＨＤＤ２３を、コンピュータ読み取り可能な記録媒体ということがある。

通信インタフェース２７は、ネットワーク２７ａに接続され、ネットワーク２７ａを介して他の情報処理装置と通信を行うインタフェースである。通信インタフェース２７は、スイッチなどの通信装置とケーブルで接続される有線通信インタフェースでもよいし、基地局と無線リンクで接続される無線通信インタフェースでもよい。

次に、最適化装置２０の機能及び処理手順を説明する。
図４は、第２の実施の形態の最適化装置の機能例を示すブロック図である。
最適化装置２０は、記憶部３０、処理部３１を有する。処理部３１は、制御部３１ａ、設定読込部３１ｂ、問題設定部３１ｃ、ハミルトニアン計算部３１ｄ、エネルギー差計算部３１ｅ、乱数発生部３１ｆ、遷移確率計算部３１ｇ、ＭＣＭＣ計算実行部３１ｈ、ＳＡ実行部３１ｉ、レプリカ交換実行部３１ｊを有する。さらに、処理部３１は、エネルギー更新部３１ｋ、結果出力部３１ｌ、λアニール部３１ｍ、パラメータ最適化部３１ｎを有する。

なお、記憶部３０は、たとえば、ＨＤＤ２３に確保した記憶領域を用いて実装できる。処理部３１は、たとえば、ＣＰＵ２１が実行するプログラムモジュールを用いて実装できる。

記憶部３０は、エネルギー情報、スピン情報、レプリカ情報、温度情報、問題設定情報、ハミルトニアン情報を記憶する。
エネルギー情報は、計算されたエネルギーの初期値や、これまで計算されたエネルギーを最小値から小さい順にＮ_ｒａｎｋ個含む。また、エネルギー情報は、Ｎ_ｒａｎｋ個のエネルギーに対応した各状態変数の値のＮ_ｒａｎｋ個の組合せを含んでいてもよい。スピン情報は、各状態変数の値を含む。レプリカ情報は、レプリカ交換法を実行するために用いられる情報であり、レプリカ数などを含む。温度情報は、温度パラメータの値（以下、単に温度という）または逆温度を含む。問題設定情報は、たとえば、使用する最適化方法（ＳＡ法、レプリカ交換法など）、レプリカ交換やアニーリングの頻度、サンプリング頻度、温度変更スケジュール、使用する遷移確率分布の情報、計算終了条件となるＭＣＭＣ計算の計算回数などを含む。ハミルトニアン情報は、たとえば、式（１）に示したエネルギー関数の重み係数（Ｗ_ｉｊ）、バイアス係数（ｂ_ｉ）、定数（Ｃ）などを含む。

制御部３１ａは、処理部３１の各部を制御して最小値求解処理を行う。
設定読込部３１ｂは、記憶部３０から上記の各種情報を、制御部３１ａが理解可能な形式で読み込む。

問題設定部３１ｃは、状態変数の初期値や、エネルギーの最小値の初期値を計算する。
ハミルトニアン計算部３１ｄは、エネルギーを計算する。
エネルギー差計算部３１ｅは、状態変数の値が反転されることによるエネルギー差を計算する。

乱数発生部３１ｆは、反転対象の状態変数をランダムに選択するための乱数を生成する。乱数発生部３１ｆとしてはどのようなものを使用してもよいが、ＭＣＭＣ計算の計算回数よりも十分大きい周期をもつメルセンヌツイスタなどの擬似乱数生成器などを用いることが望ましい。

遷移確率計算部３１ｇは、選択された遷移確率分布の式と、温度または逆温度と、計算されたエネルギー差に基づいて、遷移確率分布の値を計算する。
ＭＣＭＣ計算実行部３１ｈは、計算された遷移確率分布の値にしたがって、選択された状態変数の値を更新するか否かを判定するなどのＭＣＭＣ計算を行う。

ＳＡ実行部３１ｉは、ＳＡ法を実行するために、温度変更スケジュールに基づいて、温度を下げていく。
レプリカ交換実行部３１ｊは、レプリカ交換法を実行するため、後述の交換確率に基づいて、レプリカ間で温度を交換する。

エネルギー更新部３１ｋは、ハミルトニアン計算部３１ｄが計算したエネルギーに基づいて、記憶部３０に記憶されているエネルギー情報を更新する。
結果出力部３１ｌは、ＭＣＭＣ計算の現在の計算回数が、計算終了条件となる所定の回数（以下、Ｎ_Ｅとする）に達した場合に、エネルギー情報を計算結果として出力する。なお、結果出力部３１ｌは、Ｎ_Ｅよりも少ない計算回数（以下Ｎ_Ｉとする）のたびに、その時点で計算されたエネルギーと、そのエネルギーに対応する各状態変数の値を出力してもよい。その場合、最適化装置２０は、サンプリング装置として機能し、Ｎ_Ｉはサンプリング頻度となる。

λアニール部３１ｍは、ＳＡ法とは異なるアニーリング法である後述のλアニール法（パラメータアニール法、２状態遷移法とも呼ばれる）を実行する。
パラメータ最適化部３１ｎは、後述の方法により、各種パラメータを最適化する処理を行う。

図５は、第２の実施の形態の最適化装置の制御方法の一例の処理の流れを示すフローチャートである。
処理が開始すると、まず、設定読込部３１ｂが、記憶部３０から上記の各種情報を、制御部３１ａが理解可能な形式で読み込む（ステップＳ１０）。読み込まれる情報の例については後述する。その後、問題設定部３１ｃは、状態変数の初期化とエネルギーの初期値の計算を含む初期化処理を行う（ステップＳ１１）。そして、制御部３１ａは、たとえば、記憶部３０から読み込まれた問題設定情報にパラメータ最適化処理を行うことを指示する情報が含まれているか否かを判定する（ステップＳ１２）。制御部３１ａが、パラメータ最適化処理を行うことを指示する情報が含まれていると判定した場合、ステップＳ２４の処理が行われる。

制御部３１ａが、パラメータ最適化処理を行うことが指定されていないと判定した場合、制御部３１ａは、現在のＭＣＭＣ計算の計算回数（以下ステップ数と呼ぶ場合もある）を示すＮ_Ｃを１に初期化する（ステップＳ１３）。ステップＳ１３の処理後、ＭＣＭＣ計算実行部３１ｈは、ＭＣＭＣ計算により、Ｍ－Ｈ計算処理を行う（ステップＳ１４）。Ｍ－Ｈ計算処理の例については後述する。その後、エネルギー更新部３１ｋは、Ｍ－Ｈ計算処理の結果にしたがってハミルトニアン計算部３１ｄが計算したエネルギーに基づいて、記憶部３０に記憶されているエネルギー情報を更新するエネルギー更新処理を行う（ステップＳ１５）。

ステップＳ１５の処理後、制御部３１ａは、記憶部３０から読み込まれた問題設定情報にアニーリング（ＳＡ法またはλアニール法）を行うことを指示する情報が含まれているか否かを判定する（ステップＳ１６）。制御部３１ａが、アニーリングを行うことを指示する情報が含まれていると判定した場合、ＳＡ実行部３１ｉまたはλアニール部３１ｍによるアニーリング処理が行われる（ステップＳ１７）。アニーリング処理の例については後述する。アニーリング処理後、または、制御部３１ａが、問題設定情報にアニーリングを行うことを指示する情報が含まれていないと判定した場合、ステップＳ１８の処理が行われる。

ステップＳ１８の処理では、制御部３１ａは、記憶部３０から読み込まれた問題設定情報にレプリカ交換法を行うことを指示する情報が含まれているか否かを判定する（ステップＳ１８）。制御部３１ａが、レプリカ交換法を行うことを指示する情報が含まれていると判定した場合、レプリカ交換実行部３１ｊによるレプリカ交換処理が行われる（ステップＳ１９）。レプリカ交換処理の例については後述する。レプリカ交換処理後、または、制御部３１ａが、問題設定情報にレプリカ交換法を行うことを指示する情報が含まれていないと判定した場合、ステップＳ２０の処理が行われる。なお、ＳＡ法によるアニーリング処理とレプリカ交換処理とは互いに排他処理であるため、どちらか一方の選択が可能である。一方、λアニール法によるアニーリング処理とレプリカ交換処理とは同時に選択可能である。

ステップＳ２０の処理（サンプル出力処理）では、結果出力部３１ｌは、Ｎ_Ｉ回の計算処理を行うたびに、その時点で計算されたエネルギーと、そのエネルギーに対応する各状態変数の値を出力する。サンプル出力処理の例については後述する。

次に、制御部３１ａは、Ｎ_Ｃ＝Ｎ_Ｅであるか否かを判定する（ステップＳ２１）。制御部３１ａは、Ｎ_Ｃ＝Ｎ_Ｅではないと判定した場合、Ｎ_Ｃを＋１する（ステップＳ２２）。その後、ステップＳ１４からの処理が繰り返される。制御部３１ａがＮ_Ｃ＝Ｎ_Ｅであると判定した場合、結果出力部３１ｌは、エネルギー情報を計算結果として出力する（ステップＳ２３）。その後、最適化装置２０の処理が終了する。出力されるエネルギー情報は、たとえば、これまで計算されたエネルギーを最小値から小さい順に並べたＮ_ｒａｎｋ個のエネルギーを含む。また、エネルギー情報は、Ｎ_ｒａｎｋ個のエネルギーに対応した各状態変数の値のＮ_ｒａｎｋ個の組合せを含む。

ステップＳ２４の処理では、パラメータ最適化部３１ｎは、後述の方法により、各種パラメータを最適化する処理を行い、その結果を表示する。その後、最適化装置２０の処理が終了する。

なお、図５に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
（情報読込処理の例）
図６は、情報読込処理の一例の処理の流れを示すフローチャートである。

設定読込部３１ｂは、記憶部３０から計算終了条件となる計算回数（Ｎ_Ｅ）、サンプリング頻度（Ｎ_Ｉ）を読み込む（ステップＳ３０，Ｓ３１）。さらに、設定読込部３１ｂは、記憶部３０から温度（たとえば、前述の式（１４）～式（１６）のβ（＝１／Ｔ）を決めるＴ（温度（絶対温度）））を読み込む（ステップＳ３２）。また、設定読込部３１ｂは、記憶部３０から問題サイズ（状態変数の数（Ｎ））、ハミルトニアン情報（式（１）に示したエネルギー関数の重み係数（Ｗ_ｉｊ）、バイアス係数（ｂ_ｉ）、定数（Ｃ））を読み込む（ステップＳ３３，Ｓ３４）。

また、設定読込部３１ｂは、記憶部３０からスピン初期化法（状態変数の初期値の決め方）、遷移確率分布情報を読み込む（ステップＳ３５，Ｓ３６）。遷移確率分布情報は、ＭＣＭＣ計算にどの遷移確率分布を用いるかを示す情報（分布名など）や、使用する遷移確率分布のパラメータ（たとえば、式（１４）～式（１６）のｍ_１～ｍ_３）を含む。

さらに、設定読込部３１ｂは、記憶部３０から、使用する最適化方法（ＳＡ法、レプリカ交換法など）、アニーリング情報（たとえば、温度変更スケジュール、アニーリングの頻度）を読み込む（ステップＳ３７，Ｓ３８）。また、設定読込部３１ｂは、記憶部３０から、レプリカ数、レプリカ交換用の温度情報（各レプリカに設定する温度列）、アニーリングやレプリカ交換の頻度を読み込み（ステップＳ３９，Ｓ４０，Ｓ４１）、情報読込処理を終える。

なお、図６に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
（初期化処理の例）
イジングモデルにおいて、スピン（状態変数）の初期化は重要である。以下に示す初期化処理の例では、問題設定部３１ｃは、４通りのスピン初期化法（“Ｅｘｔｅｒｎａｌ”、“Ｒａｎｄｏｍ”、“Ｏｎｅ”、“Ｚｅｒｏ”）の何れかを行うものとしている。

図７は、初期化処理の一例の処理の流れを示すフローチャートである。
たとえば、まず問題設定部３１ｃは、スピン初期化法が“Ｅｘｔｅｒｎａｌ”であるか否かを判定する（ステップＳ５０）。問題設定部３１ｃは、スピン初期化法が“Ｅｘｔｅｒｎａｌ”であると判定した場合、最適化装置２０の外部から取得した状態変数をスピン情報の初期値として設定する（ステップＳ５１）。

問題設定部３１ｃは、スピン初期化法が“Ｅｘｔｅｒｎａｌ”ではないと判定した場合、スピン初期化法が“Ｒａｎｄｏｍ”であるか否かを判定する（ステップＳ５２）。問題設定部３１ｃは、スピン初期化法が“Ｒａｎｄｏｍ”であると判定した場合、各状態変数に０または１をランダムに設定してスピン情報の初期値とする（ステップＳ５３）。たとえば、問題設定部３１ｃは、全ての状態変数を確率５０％で０に設定し、０に設定されなかった状態変数を１に設定する。

問題設定部３１ｃは、スピン初期化法が“Ｒａｎｄｏｍ”ではないと判定した場合、スピン初期化法が“Ｏｎｅ”であるか否かを判定する（ステップＳ５４）。問題設定部３１ｃは、スピン初期化法が“Ｏｎｅ”であると判定した場合、全状態変数を１に設定してスピン情報の初期値とする（ステップＳ５５）。

問題設定部３１ｃは、スピン初期化法が“Ｏｎｅ”ではないと判定した場合、スピン初期化法が“Ｚｅｒｏ”であると判定し、全状態変数を０に設定してスピン情報の初期値とする（ステップＳ５６）。

ステップＳ５１，Ｓ５３，Ｓ５５，Ｓ５６の処理後、問題設定部３１ｃは、エネルギー（Ｅ）の初期値を計算する（ステップＳ５７）。問題設定部３１ｃは、エネルギーの初期値を、ハミルトニアン情報と、スピン情報の初期値に基づいて、式（１）により計算する。

さらに、問題設定部３１ｃは、Ｅ_ｍｉｎ＝Ｅとすることで、計算したエネルギーの初期値をエネルギーの最小値（Ｅ_ｍｉｎ）として設定し（ステップＳ５８）、初期化処理を終える。

なお、図７に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
（Ｍ－Ｈ計算処理の例）
図８は、Ｍ－Ｈ計算処理の一例の処理の流れを示すフローチャートである。

乱数発生部３１ｆは、Ｎ個の状態変数の識別番号を１～Ｎとしたとき、１～Ｎの整数値をランダムに発生することで、反転候補の状態変数を決定する（ステップＳ６０）。
そして、エネルギー差計算部３１ｅは、反転候補の状態変数の値を反転させる場合のエネルギー差（ΔＥ）を計算する（ステップＳ６１）。反転候補の状態変数の識別番号がｋである場合、エネルギー差は、ΔＥ_ｋ＝－（１－２ｘ_ｋ）ｈ_ｋである。ｈ_ｋは、式（２）により表されるローカルフィールドである。

ＭＣＭＣ計算実行部３１ｈは、エネルギー差計算部３１ｅが計算したエネルギー差が０より小さいか否かを判定する（ステップＳ６２）。ＭＣＭＣ計算実行部３１ｈがエネルギー差は０以上であると判定した場合、ＭＣＭＣ計算実行部３１ｈ（または乱数発生部３１ｆ）は、区間［０，１］で一様乱数Ｒを発生する（ステップＳ６３）。そして、ＭＣＭＣ計算実行部３１ｈは、遷移確率計算部３１ｇが計算した遷移確率分布（ｆ（ΔＥ））が、一様乱数Ｒ以上であるか否かを判定する（ステップＳ６４）。

ＭＣＭＣ計算実行部３１ｈは、ｆ（ΔＥ）≧Ｒであると判定した場合、またはステップＳ６２の処理でΔＥ＜０であると判定した場合、反転候補の状態変数の値を反転させる（ステップＳ６５）。このようにして記憶部３０に記憶されるスピン情報が更新される。

また、ハミルトニアン計算部３１ｄは、ステップＳ６１の処理で計算されたΔＥを用いて、エネルギーを計算（更新）する（ステップＳ６６）。
ＭＣＭＣ計算実行部３１ｈは、ステップＳ６４の処理で、ｆ（ΔＥ）＜Ｒであると判定した場合、またはステップＳ６６の処理後、１回のＭ－Ｈ計算処理を終える。

なお、図８に示した処理の流れは一例である。たとえば、ステップＳ６５，Ｓ６６の処理の順序を入れ替えてもよい。
（エネルギー更新処理の一例）
図９は、エネルギー更新処理の一例の処理の流れを示すフローチャートである。

エネルギー更新部３１ｋは、ハミルトニアン計算部３１ｄが計算したエネルギー（Ｅ）が、エネルギー情報として記憶部３０に記憶されているＮ_ｒａｎｋ個のエネルギーのうち、何番目に小さいのかを示すｒａｎｋをＮ_ｒａｎｋ＋１に初期化する（ステップＳ７０）。

そして、エネルギー更新部３１ｋは、Ｎ_ｒａｎｋ個のエネルギーの大きさの順位を示すｊを１とした後（ステップＳ７１）、Ｅ＜Ｅ_ｍｉｎ［ｊ］であるか否かを判定する（ステップＳ７２）。Ｅ_ｍｉｎ［ｊ］は、Ｎ_ｒａｎｋ個のエネルギーのうち、ｊ番目に小さいエネルギーを表す。

エネルギー更新部３１ｋは、Ｅ＜Ｅ_ｍｉｎ［ｊ］であると判定した場合、ｒａｎｋ＝ｊとする（ステップＳ７３）。つまり、今回計算されたエネルギーが、Ｎ_ｒａｎｋ個のエネルギーのうちに入る（ランキング内である）場合、その順位がｒａｎｋ＝ｊとして、たとえば、ＲＡＭ２２に保持される。

エネルギー更新部３１ｋは、Ｅ≧Ｅ_ｍｉｎ［ｊ］であると判定した場合、ｊ＝Ｎ_ｒａｎｋであるか否かを判定し（ステップＳ７４）、ｊ＝Ｎ_ｒａｎｋではないと判定した場合、ｊ＝ｊ＋１（ステップＳ７５）とし、ステップＳ７２からの処理を繰り返す。

エネルギー更新部３１ｋは、ｒａｎｋ＝Ｎ_ｒａｎｋ＋１のまま変わらずに、ｊ＝Ｎ_ｒａｎｋと判定した場合（エネルギーがランキング外である場合）、エネルギー更新処理を終了する。

エネルギー更新部３１ｋは、ステップＳ７３の処理後、ｒａｎｋ＝Ｎ_ｒａｎｋであるか否かを判定する（ステップＳ７６）。エネルギー更新部３１ｋは、ｒａｎｋ＝Ｎ_ｒａｎｋではないと判定した場合、ｊ＝ｒａｎｋ＋１とする（ステップＳ７７）。

そして、エネルギー更新部３１ｋは、Ｅ_ｍｉｎ［ｊ］＝Ｅ_ｍｉｎ［ｊ－１］とし（ステップＳ７８）、｛ｘ_ｊ｝＝｛ｘ_ｊ－１｝とする（ステップＳ７９）。｛ｘ_ｊ｝は、ｊ番目に小さいエネルギーに対応した各状態変数の値を示す。

その後、エネルギー更新部３１ｋは、ｊ＝Ｎ_ｒａｎｋであるか否かを判定し（ステップＳ８０）、ｊ＝Ｎ_ｒａｎｋではないと判定した場合、ｊ＝ｊ＋１とし（ステップＳ８１）、ステップＳ７８からの処理を繰り返す。

エネルギー更新部３１ｋは、ｊ＝Ｎ_ｒａｎｋであると判定した場合、または、ステップＳ７６の処理で、ｒａｎｋ＝Ｎ_ｒａｎｋであると判定した場合、Ｅ_ｍｉｎ［ｒａｎｋ］＝Ｅとし（ステップＳ８２）、｛ｘ_ｒａｎｋ｝＝｛ｘ｝とする（ステップＳ８３）。このようにして記憶部３０に記憶されているエネルギー情報が更新される。ステップＳ８３の処理後、エネルギー更新処理が終了する。

なお、図９に示した処理の流れは一例である。適宜処理の順序が入れ替えられていてもよい。
（アニーリング処理の一例）
図１０は、アニーリング処理の一例の処理の流れを示すフローチャートである。

制御部３１ａは、記憶部３０から読み込まれた問題設定情報にＳＡ法を行うことを指示する情報が含まれているか否かを判定する（ステップＳ９０）。制御部３１ａが、問題設定情報にＳＡ法を行うことを指示する情報が含まれていると判定した場合、ＳＡ実行部３１ｉは、ＳＡ処理を行う（ステップＳ９１）。制御部３１ａは、問題設定情報にＳＡ法を行うことを指示する情報が含まれていないと判定した場合、問題設定情報にλアニール法を行うことを指示する情報が含まれているか否かを判定する（ステップＳ９２）。制御部３１ａが、問題設定情報にλアニール法を行うことを指示する情報が含まれていると判定した場合、λアニール部３１ｍは、λアニール処理を行う（ステップＳ９３）。

ステップＳ９１，Ｓ９３の処理後、または、ステップＳ９２の処理で問題設定情報にλアニール法を行うことを指示する情報が含まれていないと判定された場合、アニーリング処理が終了する。

（ＳＡ処理の一例）
図１１は、ＳＡ処理の一例の処理の流れを示すフローチャートである。
ＳＡ実行部３１ｉは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ａ）＝０であるか否かを判定する（ステップＳ１００）。Ｎ_Ｃは現在の計算回数（ステップ数）であり、Ｎ_Ａは、アニーリングの頻度を示すステップ数である。ステップＳ１００の処理では、Ｎ_ＣがＮ_Ａの倍数であるか否かが判定される。ＳＡ実行部３１ｉは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ａ）＝０でないと判定した場合、ＳＡ処理を終了する。

ＳＡ実行部３１ｉは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ａ）＝０であると判定した場合、温度変更スケジュールを判定する（ステップＳ１０１）。
図１１の例では、４つの温度変更スケジュールの例が示されている。“ｉｎｖ”は、温度（Ｔ）を、Ｎ_Ｃに反比例させて減少させる温度変更スケジュールである。“ｉｎｖｒｏｏｔ”は、Ｔの値を、Ｎ_Ｃの平方根に反比例させて減少させる温度変更スケジュールである。“ｅｘｐ”は、Ｔの値を、指数型スケジュールで減少させる温度変更スケジュールである。“ｌｉｓｔ”は、Ｔの値を、リスト（任意に与えられたステップ数とＴの値との組）にしたがって変化（減少）させる温度変更スケジュールである。

温度変更スケジュールが“ｉｎｖ”である場合、ＳＡ実行部３１ｉは、Ｔの値をＮ_Ｃに反比例させて減少させる（ステップＳ１０２）。温度変更スケジュールが“ｉｎｖｒｏｏｔ”である場合、ＳＡ実行部３１ｉは、Ｔの値をＮ_Ｃの平方根に反比例させて減少させる（ステップＳ１０３）。温度変更スケジュールが“ｅｘｐ”である場合、ＳＡ実行部３１ｉは、Ｔを指数型スケジュールで減少させる（ステップＳ１０４）。温度変更スケジュールが“ｌｉｓｔ”である場合、ＳＡ実行部３１ｉは、Ｔの値をリストにしたがって減少させる（ステップＳ１０２）。リストを用いることで、任意の温度変更スケジュールを設定できる。

ステップＳ１０２～Ｓ１０５の処理後、ＳＡ処理が終了する。
（λアニール処理の一例）
λアニール法では、エネルギー関数として、Ｈ（λ）＝（１－λ）Ｈ_Ａ＋λＨ_Ｂが用いられる。

Ｈ_Ｂは、答えを求めたいハミルトニアンであり、たとえば、式（１）が用いられる。Ｈ_Ａは初期状態のハミルトニアンであり、答えが良く知られた既知のハミルトニアンである。たとえば、１次元的なスピンの並びを考え、隣接するスピン間のみ相互作用をし、相互作用の強さは全て一定の値として与えられるものとする。こうすると、Ｈ_Ａには解析解が存在する状態を作り出せる。たとえば、Ｈ_Ａとして、全てのスピン（状態変数）の値が０のとき、エネルギー（Ｈ（λ））が最低の状態になるようなハミルトニアンを採用することができる。もちろん、全ての状態変数の値が１のとき、エネルギーが最低の状態になるようなハミルトニアンをＨ_Ａとして採用してもよい。Ｈ_Ａの選び方は、初期状態の決め方に依存する。

λ（アニーリング変数）は、０≦λ≦１であり、時間の関数としてλ＝λ（ｔ）として定義することができる。ただし、λ（ｔ）は、計算の開始時刻（ｔ＝０）にλ（０）＝０、計算の終了時刻（ｔ＝τ）にλ（τ）＝１という境界条件を満たすものとする。このため、時間の関数となるハミルトニアンは、計算の開始時刻では、Ｈ（０）＝Ｈ_Ａとなり、計算の終了時刻では、Ｈ（τ）＝Ｈ_Ｂとなる。

したがって、ｔ＝０にて全ての状態変数が初期化された状態から計算が開始され、計算終了時に答えの知りたいハミルトニアンの答えが得られる。これはイジング型量子コンピュータの考え方であるが、最適化装置２０は、量子論を扱わず、古典的なハミルトニアンの範囲内でエネルギーの最小値を求める手法として、λアニール法を用いる。

図１２は、λアニール処理の一例の処理の流れを示すフローチャートである。
λアニール部３１ｍは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ａ）＝０であるか否かを判定する（ステップＳ１１０）。Ｎ_Ｃは現在の計算回数（ステップ数）であり、Ｎ_Ａは、アニーリングの頻度を示すステップ数である。ステップＳ１１０の処理では、Ｎ_ＣがＮ_Ａの倍数であるか否かが判定されている。λアニール部３１ｍは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ａ）＝０でないと判定した場合、λアニール処理を終了する。

λアニール部３１ｍは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ａ）＝０であると判定した場合、たとえば、読み込まれた問題設定情報に含まれるλ変更スケジュールを判定する（ステップＳ１１１）。

図１２の例では、４つのλ変更スケジュールの例が示されている。“ｌｉｎｅａｒ”は、λ（Ｎ_Ｃ）＝Ｎ_Ｃ／Ｎ_Ｅにしたがって、λ（Ｎ_Ｃ）をＮ_Ｃの増加とともに増加させるλ変更スケジュールである。“ｑｕａｄｒａｔｉｃ”は、λ（Ｎ_Ｃ）を（Ｎ_Ｃ／Ｎ_Ｅ）の２乗にしたがって、Ｎ_Ｃの増加とともに増加させるλ変更スケジュールである。“ｓｑｒｔ”は、λ（Ｎ_Ｃ）を（Ｎ_Ｃ／Ｎ_Ｅ）の平方根にしたがって、Ｎ_Ｃの増加とともに増加させるλ変更スケジュールである。“ｌｉｓｔ”は、λ（Ｎ_Ｃ）をｆ（Ｎ_Ｃ）という関数にしたがって増加させるλ変更スケジュールである。ｆ（Ｎ_Ｃ）は任意であるが、ステップ数とそのステップ数におけるλの値による組（Ｎ_ｉ，λ_ｉ）の列として与えられる。

λ変更スケジュールが“ｌｉｎｅａｒ”である場合、λアニール部３１ｍは、λ（Ｎ_Ｃ）を、λ（Ｎ_Ｃ）＝Ｎ_Ｃ／Ｎ_Ｅを計算することで求める（ステップＳ１１２）。λ変更スケジュールが“ｑｕａｄｒａｔｉｃ”である場合、λアニール部３１ｍは、λ（Ｎ_Ｃ）を、（Ｎ_Ｃ／Ｎ_Ｅ）の２乗を計算することで求める（ステップＳ１１３）。λ変更スケジュールが“ｓｑｒｔ”である場合、λアニール部３１ｍは、λ（Ｎ_Ｃ）を、（Ｎ_Ｃ／Ｎ_Ｅ）の平方根を計算することで求める（ステップＳ１１４）。λ変更スケジュールが“ｌｉｓｔ”である場合、λアニール部３１ｍは、λ（Ｎ_Ｃ）を、前述のｆ（Ｎ_Ｃ）という関数にしたがって計算する（ステップＳ１１５）。

ステップＳ１１２～Ｓ１１５の処理後、λアニール部３１ｍは、変更前のλについてのＨ（λ）（たとえば、ハミルトニアン計算部３１ｄによって計算される）の最小値を与える各状態変数の値を記憶部３０から読み込む（ステップＳ１１６）。読み込まれた各状態変数の値は、λ変更後に再度Ｍ－Ｈ計算処理が行われるときの初期値として用いられる。ステップＳ１１６の処理後、λアニール処理が終了する。

なお、λアニール処理は、以下に示すレプリカ交換処理と同時に実行させることができる。
（レプリカ交換処理の一例）
図１３は、レプリカ交換処理の一例の処理の流れを示すフローチャートである。

レプリカ交換実行部３１ｊは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ｒ）＝０であるか否かを判定する（ステップＳ１２０）。Ｎ_Ｃは現在の計算回数（ステップ数）であり、Ｎ_Ｒは、レプリカ交換の頻度を示すステップ数である。ステップＳ１２０の処理では、Ｎ_ＣがＮ_Ｒの倍数であるか否かが判定されている。

レプリカ交換実行部３１ｊは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ｒ）＝０であると判定した場合、ＭＯＤ（Ｎ_Ｃ，２）＝０であるか否かを判定する（ステップＳ１２１）。ステップＳ１２１の処理では、Ｎ_Ｃが偶数であるか否かが判定されている。

レプリカ交換実行部３１ｊは、ＭＯＤ（Ｎ_Ｃ，２）＝０であると判定した場合、温度順が偶数番目のレプリカと１つ温度が大きいレプリカによるペアを全て選択する（ステップＳ１２２）。レプリカ交換実行部３１ｊは、ＭＯＤ（Ｎ_Ｃ，２）＝０でないと判定した場合、温度順が奇数番目のレプリカと１つ温度が大きいレプリカによるペアを全て選択する（ステップＳ１２３）。このようにレプリカ交換候補のレプリカペアを決めることで、レプリカ交換が効率化される。

ステップＳ１２２，Ｓ１２３の処理後、レプリカ交換実行部３１ｊは、レプリカ交換を行うペアのレプリカ番号を決定する（ステップＳ１２４）。
図１４は、レプリカ交換を行うペアのレプリカ番号の決定方法の一例を示す図である。

図１４において、カッコ内の上段がレプリカ番号を示し、下段がレプリカ番号のレプリカに設定されている温度を示す。ただし、Ｔ_１＜Ｔ_２＜Ｔ_３である。レプリカ交換において交換対象のレプリカ間で温度の交換を行う場合、レプリカ番号と温度との対応関係が変化する。レプリカ交換を隣接する温度が設定されたレプリカ間で行う場合、そのようなレプリカを検出するために、ステップＳ１２４の処理では、隣接する温度（Ｔ_ｉとＴ_ｉ＋１）に対応するレプリカ番号が検出される。

図１４の、一番左側の例では、隣接する温度（Ｔ_１，Ｔ_２）に対応するレプリカ番号＝１，２が検出されている。レプリカ番号＝１，２のレプリカ間で温度が交換された場合、図１４の真ん中の例のようになり、隣接する温度に対応するレプリカ番号は、１，２と、１，３の２組がある。レプリカ番号＝１，２の間では同じ温度が交換されたばかりであるので、レプリカ番号＝１，３がレプリカ交換を行うペアのレプリカ番号として決定される。

なお、レプリカ間で状態変数の値を交換する場合には、このような処理が不要であるが、状態変数の自由度が高くなると交換時に送受信より発生する情報量が多くなるため、上記のように温度を交換することが望ましい。

その後、レプリカ交換実行部３１ｊは、レプリカ交換を行うペアのｉｄ（ｊ＝０～Ｎ_ｒｅｐ－１（Ｎ_ｒｅｐはステップＳ１２２，Ｓ１２３で選択されたペア数））を０にする（ステップＳ１２５）。そして、レプリカ交換実行部３１ｊは、ｉｄ＝ｊのペアに対して、レプリカ交換を実行するか否かを判定する（ステップＳ１２６）。

レプリカ交換法における交換確率は、用いる遷移確率分布がボルツマン分布のように詳細つり合いの原理を満たす場合、以下の式（１８）で表せる。

式（１８）において、Ｐ_Ａ（ｔ）は、レプリカ交換前の状態Ａが実現する確率である。状態Ａは、β_ｉが設定されているレプリカ番号＝ｉのレプリカにおける状態変数が確率分布π（β_ｉ，｛ｘ_ｉ｝）にしたがい、β_ｊが設定されているレプリカ番号＝ｊのレプリカにおける状態変数が確率分布π（β_ｊ，｛ｘ_ｊ｝）にしたがう状態である。β_ｉ，β_ｊは逆温度である。式（１８）において、Ｐ_Ｂ（ｔ）は、レプリカ交換後の状態Ｂが実現する確率である。状態Ｂは、β_ｉが設定されているレプリカ番号＝ｉのレプリカにおける状態変数が確率分布π（β_ｉ，｛ｘ_ｊ｝）にしたがい、β_ｊが設定されているレプリカ番号＝ｊのレプリカにおける状態変数が確率分布π（β_ｊ，｛ｘ_ｉ｝）にしたがう状態である。式（１８）では、状態変数｛ｘ_ｉ｝、｛ｘ_ｊ｝がレプリカ番号＝ｉ，ｊのレプリカ間で交換される場合が示されているが、β_ｉ，β_ｊが交換される場合も同じ交換確率Ｐ_Ａ→Ｂとなる。

一方、用いる遷移確率分布が詳細つり合いを満たさない場合、交換確率Ｐ_Ａ→Ｂは、以下の式（１９）で表される。

式（１９）において、ΔＥはＥ_ｊ－Ｅ_ｉである。すなわち、ΔＥは、β_ｊが設定されているレプリカ番号＝ｊのレプリカにおけるエネルギーと、β_ｉが設定されているレプリカ番号＝ｉのレプリカにおけるエネルギーとの差分である。

式（１９）を用いる理由を説明する。前述の式（９）により遷移確率を定義した場合、前述のように、固有値１をもつ確率行列が得られ、定常マルコフ連鎖が形成される。つまり、状態は定常状態に収束する。この定常状態に収束したときの状態密度は、遷移確率に対応した状態密度になる。いま、レプリカ交換時に全てのレプリカ同士が、式（１９）の交換確率で交換される場合、全レプリカによる系にも定常状態が存在する。レプリカ交換法などの拡張アンサンブル法を用いる場合、マルコフ連鎖が定常状態を形成することが前提であるため、式（１９）のような定常状態が得られる交換確率が用いられる。

式（１９）のｆ（β_ｉΔＥ）として、式（１４）～式（１６）に示したような遷移確率分布が用いられる。
ステップＳ１２６の処理では、レプリカ交換実行部３１ｊは、ｉｄ＝ｊのペアに対して、式（１９）の交換確率にしたがって、レプリカ交換を実行するか否かを判定する。たとえば、レプリカ交換実行部３１ｊは、式（１９）のｆ（β_ｉΔＥ）と、０≦Ｒ≦１の一様乱数Ｒとの比較結果に基づいて、ｆ（β_ｉΔＥ）≧Ｒであるならば、レプリカ交換実行部３１ｊは、レプリカ交換を実行すると判定する。ｆ（β_ｉΔＥ）≧Ｒでないならば、レプリカ交換実行部３１ｊは、レプリカ交換を実行しないと判定する。

レプリカ交換実行部３１ｊは、レプリカ交換を実行すると判定した場合、ｉｄ＝ｊのペアに含まれるレプリカ間で温度を交換することで、レプリカ交換を実行する（ステップＳ１２７）。

レプリカ交換実行部３１ｊは、レプリカ交換を実行しないと判定した場合、またはステップＳ１２７の処理後、ｊ＝Ｎ_ｒｅｐ－１であるか否かを判定する（ステップＳ１２８）。レプリカ交換実行部３１ｊは、ｊ＝Ｎ_ｒｅｐ－１であると判定した場合、レプリカ交換処理を終了し、ｊ＝Ｎ_ｒｅｐ－１ではないと判定した場合、ｊ＝ｊ＋１とし（ステップＳ１２９）、ステップＳ１２６からの処理を繰り返す。

（サンプル出力処理の一例）
図１５は、サンプル出力処理の一例の処理の流れを示すフローチャートである。図１５では、レプリカ交換処理が行われる場合の、サンプル出力処理の例が示されている。

結果出力部３１ｌは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ｉ）＝０であるか否かを判定する（ステップＳ１３０）。Ｎ_Ｃは現在の計算回数（ステップ数）であり、Ｎ_Ｉは、サンプリング頻度を示すステップ数である。結果出力部３１ｌは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ｉ）＝０ではないと判定した場合、サンプリング出力処理を終了する。

結果出力部３１ｌは、ＭＯＤ（Ｎ_Ｃ，Ｎ_Ｉ）＝０ではあると判定した場合、レプリカ番号（ｒ）を１に設定する（ステップＳ１３１）。ステップＳ１３１の処理後、結果出力部３１ｌは、レプリカ番号＝ｒのレプリカの現在のエネルギー（Ｅ_ｒ）を出力し（ステップＳ１３２）、そのレプリカの現在の各状態変数の値（｛ｘ_ｒ｝）を出力する（ステップＳ１３３）。なお、結果出力部３１ｌは、λアニール法が実行されている場合、エネルギーとして各λの値についてのＨ（λ）を出力する。

さらに、結果出力部３１ｌは、Ｎ_Ｃと、レプリカ番号＝ｒのレプリカに設定されている温度（Ｔ_ｒ）と、レプリカ番号（ｒ）とを出力する（ステップＳ１３４）。
その後、結果出力部３１ｌは、ｒ＝Ｎ_ｒであるか否かを判定する（ステップＳ１３５）。Ｎ_ｒはレプリカ数である。結果出力部３１ｌは、ｒ＝Ｎ_ｒであると判定した場合、サンプリング出力処理を終了し、ｒ＝Ｎ_ｒではないと判定した場合、ｒ＝ｒ＋１とし（ステップＳ１３６）、ステップＳ１３２からの処理を繰り返す。

（パラメータ最適化処理）
エネルギーの最小値をヒューリスティックに求める場合、計算者は最適な計算結果を最小計算回数または最小時間で知りたい。しかし、用いる遷移確率分布によっては解が極小値にトラップされ、簡単には脱出できない場合もある。その場合、結果としてエネルギー空間上を広くサンプリングできず（サンプリング効率が悪化し）、最小値または次善の解の求解精度と速度が悪化する。このため、最小値求解問題を計算する際には、遷移確率分布の選択が重要である。

図１６は、２種類の遷移確率分布を用いた場合のサンプリング結果の比較例を示す図である。図１６の左では、式（１６）で表される遷移確率分布（以下Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布と呼ぶ場合もある）であり、ｍ_３＝１／４とした場合の遷移確率分布を用いた場合のサンプリング結果が示されている。図１６の右では、式（１４）で表される遷移確率分布（以下Ｐｏｗｅｒ ｌａｗ型の遷移確率分布と呼ぶ場合もある）であり、ｍ_１≒１とした場合の遷移確率分布を用いた場合のサンプリング結果が示されている。ただし、両分布とも、温度（絶対温度）を０．１としている。

横軸はＭｏｎｔｅ Ｃａｒｌｏ Ｓｔｅｐｓ（前述のステップ数に相当する）を表し、縦軸はＣｏｓｔ Ｆｕｎｃｔｉｏｎ（前述のエネルギーに相当する）を表す。
図１６の例では、サンプリング結果の分散が大きく（サンプリング範囲が広く）、よりエネルギーが低い領域を探索できるのは、Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合である。

しかしながら、温度などのパラメータを変えると、結論が変わりえる。
図１７は、温度を変更したときの２種類の遷移確率分布を用いた場合のサンプリング結果の比較例を示す図である。図１７では、温度を１．０と変えた以外は、図１６と同じ条件のＴｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布と、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合のサンプリング結果が示されている。

温度を１．０とした場合、図１７に示されているように、よりエネルギーが低い領域を探索できるのは、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合である。一方でサンプリング結果の分散が大きいのはＴｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合である。

このように、最小値求解問題を計算する際には、最適な遷移確率分布を選択することや、最適なパラメータを推定することが重要である。
一方で、最適な遷移確率分布の選択や、最適なパラメータの推定を行うための手続きは煩雑である。最小値求解問題の本計算をする前に、使用する遷移確率分布やパラメータを変えながら多数回の計算を行うことになるためである。

そこで、パラメータ最適化部３１ｎは、遷移確率分布に含まれる最適なパラメータを、以下のように自動推定する機能を有する。なお、自動推定機能は、以下の性質を満たすことが望ましい。すなわち、短時間の予備計算で推定が行えること、サンプリング時のエネルギーの平均値が小さくなること、サンプリング時のエネルギーの分散が大きくなること、エネルギーの最小値がより小さくなること、という４つの性質である。

図１８は、パラメータ最適化処理の一例の処理の流れを示すフローチャートである。
パラメータ最適化部３１ｎは、最適化対象の遷移確率分布の関数名を、記憶部３０から読み込む（ステップＳ１４０）。たとえば、式（１４）で表される遷移確率分布の関数名は“Ｐｏｗｅｒ ｌａｗ”、式（１５）で表される遷移確率分布の関数名は“Ｈｙｐｅｒｂｏｌｉｃ”、式（１６）で表される遷移確率分布の関数名は“Ｔｅｍｐｅｒｅｄ－Ｇａｕｓｓｉａｎ”などと指定される。また、ユーザが指定する遷移確率分布が読み込まれる場合、“Ｕｓｅｒ－Ｄｅｆｉｎｅｄ”などと関数名が指定される。パラメータ最適化部３１ｎは、“Ｕｓｅｒ－Ｄｅｆｉｎｅｄ”を読み込む場合、（ΔＥ_ｉ，ｆ（ΔＥ_ｉ））の組合せも読み込み、０≦ｆ（ΔＥ_ｉ）≦１を満たすように、各組合せによる各点を滑らかに繋いでいく。このときパラメータ最適化部３１ｎは、たとえば、キュービックスプラインなどを用いることができる。なお、パラメータ最適化部３１ｎは、ΔＥ_ｉの最大値が与えられていて、読み込まれたΔＥ_ｉがその最大値より大きい場合は、ｆ（ΔＥ_ｉ）＝０とする。

そして、パラメータ最適化部３１ｎは、最適化対象のパラメータ名を、記憶部３０から読み込む（ステップＳ１４１）。全ての遷移確率分布の関数において共通であるパラメータは、温度（式（１４）～式（１６）では逆温度（β＝１／Ｔ）で表されている）であるが、式（１４）～式（１６）の場合、ｍ_１～ｍ_３も最適化対象のパラメータとして選択できる。

さらに、パラメータ最適化部３１ｎは、パラメータを最適化する範囲を、記憶部３０から読み込む（ステップＳ１４２）。たとえば、パラメータが温度の場合、最小温度と最大温度が読み込まれる。なお、記憶部３０に記憶されている範囲とは異なる範囲を指定する場合には、たとえば、入力デバイス２５ａから指定すべき範囲が入力される。たとえば、パラメータが温度の場合、本実施の形態の最適化装置２０では負の絶対温度は定義されないため、範囲は０以上となる。

次に、パラメータ最適化部３１ｎは、パラメータを変化させる回数（予備計算回数）Ｎ_ｐを、記憶部３０から読み込む（ステップＳ１４３）。そして、パラメータ最適化部３１ｎは、パラメータの刻み幅ΔＰを計算する（ステップＳ１４４）。パラメータの最大値をＰ_ｍａｘ、パラメータの最小値をＰ_ｍｉｎとしたとき、パラメータ最適化部３１ｎは、たとえば、ΔＰ＝（Ｐ_ｍａｘ－Ｐ_ｍｉｎ）／（Ｎ_ｐ－１）を計算することで、ΔＰを求める。なお、パラメータ最適化部３１ｎは、Ｐ_ｍｉｎとＰ_ｍａｘの間を不等分に分割したΔＰを示すリストを読み込んでもよい。

その後、パラメータ最適化部３１ｎは、切捨て対象の計算回数Ｎ_ｄｕｍｐを、記憶部３０から読み込む（ステップＳ１４５）。Ｎ_ｄｕｍｐは、数値統計の分野で“バーンイン”と呼ばれる期間を決める計算回数であり、マルコフ連鎖が定常状態に収束するまでのステップ数を意味する。計算対象の問題によってバーンインのステップ数は異なる。

なお、ここまでの説明では、パラメータ最適化部３１ｎは、記憶部３０から各種情報を読み込むものとしたが、図３に示したディスプレイ２４ａに読み込み候補を表示させ、ユーザによる入力デバイス２５ａの操作によって入力された情報を読み込んでもよい。

ステップＳ１４５の処理後、パラメータ最適化部３１ｎは、現在の予備計算回数を示すｉを１とし（ステップＳ１４６）、パラメータの値を計算する（ステップＳ１４７）。ステップＳ１４７の処理では、パラメータ最適化部３１ｎは、Ｐ_ｉ＝Ｐ_ｍｉｎ＋（ｉ－１）を計算することで、パラメータの値（Ｐ_ｉ）を求める。

その後、パラメータ最適化部３１ｎは、ＭＣＭＣ計算の現在の計算回数を示すｊを１とし（ステップＳ１４８）、ＭＣＭＣ計算によりエネルギー（Ｅ_ｉ）とエネルギー差（ΔＥ_ｉ）を計算する（ステップＳ１４９）。

そして、パラメータ最適化部３１ｎは、ｊ＞Ｎ_ｄｕｍｐであるか否か（計算回数がバーンイン期間を経過したか否か）を判定する（ステップＳ１５０）。パラメータ最適化部３１ｎは、ｊ＞Ｎ_ｄｕｍｐであると判定した場合、Ｅ_ｉとΔＥ_ｉとエネルギーの最小値を、たとえば、図３に示したＲＡＭ２３に保存する（ステップＳ１５１）。

パラメータ最適化部３１ｎは、ｊ＞Ｎ_ｄｕｍｐではないと判定した場合（マルコフ連鎖が定常状態に達していないと判定した場合）、またはステップＳ１５１の処理後、ｊ＝Ｍであるか否かを判定する（ステップＳ１５２）。Ｍは、ＭＣＭＣ計算の計算回数である。Ｍを大きくとりすぎると予備計算に時間がかかるため、バーンイン期間に相当する計算回数よりも大きく、バーンイン期間が経過するまでのデータを切り捨てた後のデータを用いて統計量を計算するためにかかる回数が指定される。モンテカルロ積分を収束するための計算量は１／（Ｍの平方根）に比例するため、推定を行える程度の計算回数が指定される。たとえば、Ｍは、１０００回から数百万回程度であり、記憶部３０から読み込まれる。これにより、以下に示す評価値の計算が、エネルギー関数の最小値の探索を行う際の計算時間（本計算時間）よりも短い時間で計算されるようにしている。

パラメータ最適化部３１ｎは、ｊ＝Ｍではないと判定した場合、ｊ＝ｊ＋１とし（ステップＳ１５３）、ステップＳ１４９からの処理を繰り返す。パラメータ最適化部３１ｎは、ｊ＝Ｍであると判定した場合、保存したＥ_ｉとΔＥ_ｉの平均値と標準偏差を計算する（ステップＳ１５４）。なお、計算した平均値と標準偏差は、エネルギーの最小値とともに、各パラメータの値におけるエネルギー関数の値の評価値として、たとえばＲＡＭ２３に保存される。

その後、パラメータ最適化部３１ｎは、ｉ＝Ｎ_ｐであるか否かを判定する（ステップＳ１５５）。パラメータ最適化部３１ｎは、ｉ＝Ｎ_ｐではないと判定した場合、ｉ＝ｉ＋１とし（ステップＳ１５６）、ステップＳ１４７からの処理を繰り返す。

パラメータ最適化部３１ｎは、ｉ＝Ｎ_ｐであると判定した場合、保存していた評価値を出力し（ステップＳ１５７）、パラメータ最適化処理を終える。評価値は、たとえば、図３に示したディスプレイ２４ａに出力（表示）される。なお、評価値は、記憶部３０に出力されるようにしてもよい。

なお、図１８に示した処理の流れは一例である。適宜処理の順序が入れ替えられていてもよい。
図１９は、評価値の出力例を示す図である。

図１９では、温度についての評価値のディスプレイ２４ａへの表示例が示されている。Ｔ＝１～１０までＮ_ｐ＝１０で計算を行った例が示されている。各温度についてのエネルギーの最小値（Ｅ_ｍｉｎ）、エネルギーの平均値（Ｅ_ａｖｅ）、エネルギーの標準偏差（Ｅ_{ｓｔｄｄｅｖ}）、エネルギー差の平均値（ｄＥ_ａｖｅ）、エネルギー差の標準偏差（ｄＥ_{ｓｔｄｄｅｖ}）が評価値としてリスト化されている。また、図１９の例では、Ｅ_ｍｉｎが小さい順に、各温度が順位付けされている。Ｎ_ｐが大きくなると、計算結果のデータ量が増えるため、上記のように各パラメータの値に順位付けを行うことで、よいパラメータの値の判別が容易になり、解析時間を削減できる。なお、順位付けは、Ｅ_ｍｉｎが小さい順に行うことに限定されず、Ｅ_ａｖｅが小さい順に順位付けされるようにしてもよい。

以上のような、パラメータ最適化処理を行うことで、計算者によるパラメータを決定する手間が削減され、よい計算結果が得られることが期待できるパラメータを容易に決定できる。

（効果）
以上のような第２の実施の形態の最適化装置２０によれば、第１の実施の形態の最適化装置１０と同様に以下の効果が得られる。

すなわち、ΔＥが（正に）大きいときの遷移確率がボルツマン分布よりも大きくなる遷移確率分布を適用するため、解が局所解から高効率で脱出可能となる。また、ダイナミックオフセット法のようにエネルギーにオフセットを加える方法ではないためマルコフ連鎖を破壊することもない。以上のことから、マルコフ連鎖を破壊せずに効率的にエネルギー関数の最小値の探索が可能となる。

さらに、最適化装置２０は、レプリカ交換処理の際に、上記のような遷移確率分布を用いた交換確率を採用することで、ボルツマン分布を用いる場合よりもレプリカ交換頻度が上がる。このため、少ないレプリカ数でも温度空間を各レプリカがランダムウォークするようになる。これによりサンプリング空間が広がり、サンプリング効率が向上し、より精度の高い解を得ることができる。

また、最適化装置２０は、図９に示したエネルギー更新処理において、これまで計算されたエネルギーを最小値から小さい順にＮ_ｒａｎｋ個含むようにエネルギー情報を更新していく。これにより、エネルギー関数の真の最小値が得られない場合でも、次善の最小値（解）が得られるようになる。

（計算例）
以下、上記効果を示すために、最適化装置２０による巡回セールスマン問題の計算例を示す。

巡回セールスマン問題とは、Ｎ個の都市が点Ｐ_ｉ（ｉ＝１，…，Ｎ）として定義されているとき、セールスマンが各都市を必ず１回訪問し、最後に始めの都市に戻ってくる行程の最小距離を求める問題である。ここで、点Ｐ_ｉと点Ｐ_ｊ間の距離ｄ_ｉｊは、都市間の距離を表す。巡回セールスマン問題をイジングモデルに変換するために、ｉ番目に都市ａを訪問するか否かを示す状態変数をｎ_ｉ，ａとすると、Ｎ^２個の状態変数をもつイジング型のエネルギー関数（Ｅ）は以下の式（２０）のように表せる。

しかし、式（２０）のままであると、セールスマンがなにも行動しない（全ての状態変数が０）という自明解や、セールスマンがｉ番目に全ての都市を訪問するような非現実解まで許容してしまうため、以下の式（２１）、式（２２）で表される拘束条件が課される。

式（２１）、式（２２）で表される拘束条件をペナルティ項として、式（２０）に取り込んだ以下の式（２３）が、エネルギー最小化対象の式となる。

式（２３）において、ｋ_１及びｋ_２はペナルティの大きさを表す定数である。ペナルティの大きさをエネルギーの最小値よりも大きくしておくことで、得られる自明解や非現実解はエネルギーが大きくなり、解の候補から外れる。

なお、式（２３）は、たとえば、式（１）に示した形式に変換されて計算される。ｋ_１及びｋ_２は重み係数（Ｗ_ｉｊ）に反映される。
計算対象の問題として、上記のような巡回セールスマン問題を採用したのは以下の理由による。巡回セールスマン問題は、物理の問題ではないため、ボルツマン分布に縛られる必要がないという点が１つ目の理由である。また、巡回セールスマン問題は、イジングモデルに変換できるという点が２つ目の理由である。また、式（２３）において２次の項の相互作用の強さがバラバラであり、遠く離れたスピン間（状態変数間）の相互作用も強いという点が３つ目の理由である。さらに、拘束条件がペナルティ項として式（２３）に含まれるため、ペナルティ付のイジングモデル（外場ありのイジングモデル）として理想的であるという点が４つ目の理由である。

図２０は、３種の遷移確率分布を使用した場合の定常状態の計算結果の一例を示す図である。横軸は、ステップ数を表し、縦軸はエネルギー（Ｅ）を表す。図２０では、遷移確率分布として、左から、ボルツマン分布、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布、Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合のサンプリング結果が示されている。

使用した巡回セールスマン問題は、ＴＳＰＬＩＢのｕｌｙｓｅｅｓ１６（TSPLIB, [online], ［平成３１年４月１１日検索］、インターネット<URL: https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/>参照）である。なお、計算条件は、Ｔ＝１００．０、式（２３）のｋ_１，ｋ_２は１５０、Ｎ_Ｅ＝１００，０００、Ｎ_Ｉ＝１０、式（１４）のｍ_１は１．０２、式（１６）のｍ_３は０．２５である。

図２０に示されているように、Ｐｏｗｅｒ ｌａｗ型、Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合でも、ボルツマン分布を用いた場合と同様に、マルコフ連鎖を繰り返すことによって定常状態に収束していることが分かる。さらに、エネルギーの探索範囲は、ボルツマン分布を用いた場合、［１０００，４０００］程度であるが、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布では［１０００，９０００］程度まで広がる。Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合も同様である。

図２１は、探索性能の比較例を示す図である。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。図２１では、遷移確率分布としてボルツマン分布を用いた場合のサンプリング結果（左上）、遷移確率分布としてボルツマン分布を用い、エネルギーにオフセットを加える手法（ダイナミックオフセット法）によるサンプリング結果（右上）が示されている。さらに、遷移確率分布として、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合のサンプリング結果（左下）と、Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合のサンプリング結果（右下）が示されている。

使用した巡回セールスマン問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、計算条件は、Ｔ＝０．２５、式（２３）のｋ_１，ｋ_２は１５０、Ｎ_Ｅ＝２００，０００、Ｎ_Ｉ＝１０、式（１４）のｍ_１＝は１．０２、式（１６）のｍ_３は０．２５である。なお、ダイナミックオフセット法の計算では、５０回連続して同じ状態に留まり続けた場合に、解が極小値に陥ったものとみなし、オフセットを加えるものとしている。

図２１から明らかなように、オフセットを用いずに、ボルツマン分布を遷移確率分布とした場合、一度、解が極小値に陥った場合、解が極小値から抜け出すことはない。ダイナミックオフセット法を適用した場合、解は極小値から容易に抜け出すことができるが、解として知りたいエネルギー領域のサンプリングは数点しかない（良解が少ない）。一方、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布や、Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合、極小値から抜け出すことができるとともに、より低いエネルギー領域にも多くの解が探索できている。

このため、Ｐｏｗｅｒ ｌａｗ型やＴｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合、ボルツマン分布を用いた場合やボルツマン分布を用いたダイナミックオフセット法を適用した場合よりも、探索性能が向上することが分かった。

さらに、Ｐｏｗｅｒ ｌａｗ型や、Ｔｅｍｐｅｒｅｄ Ｇａｕｓｓｉａｎ型の遷移確率分布を用いた場合は、エネルギーにオフセットを加える手法ではないためマルコフ連鎖を壊さず、レプリカ交換法のような拡張アンサンブル法をそのまま適用できる。

以下、レプリカ交換法を、遷移確率分布としてボルツマン分布を用いて実施した場合と、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いて実施した場合の計算例を示す。使用した巡回セールスマン問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、計算条件は、式（２３）のｋ_１，ｋ_２は１５０、Ｎ_Ｅ＝２００，０００（レプリカごと）、Ｎ_Ｉ＝１０、式（１４）のｍ_１＝は１．０２とした。なお、レプリカ数は８つ、各レプリカに設定される温度列は、Ｔ＝１．０，３．０，５．０，１０．０，１５．０，２０．０，２５．０，３０．０であるものとした。

図２２は、ボルツマン分布を用いた場合の各レプリカにおける温度の遷移の計算結果の一例を示す図である。図２２には、レプリカ番号＝１～８までの８つのレプリカのそれぞれにおける、温度（Ｔ）の遷移の計算結果が示されている。横軸はステップ数を表し、縦軸はＴであり、正の実数を表す。

図２２に示されているように、ボルツマン分布を用いた場合、各レプリカは温度空間を等確率でランダムウォークしているとは言えないことが分かる。たとえば、レプリカ番号＝５のレプリカは、ステップ数が２００付近から、Ｔ＝１．０に固定されてしまっている。

図２３は、ボルツマン分布を用いた場合の各レプリカにおけるエネルギーの計算結果の一例を示す図である。図２３には、レプリカ番号＝１～８までの８つのレプリカのそれぞれにおける、エネルギー（Ｅ）の計算結果が示されている。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。

図２３に示されているように、ボルツマン分布を用いた場合、ほとんどのレプリカのエネルギーは極小値に拘束され、脱出できないことが分かる。
図２４は、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合の各レプリカにおける温度の遷移の計算結果の一例を示す図である。図２４には、レプリカ番号＝１～８までの８つのレプリカのそれぞれにおける、温度（Ｔ）の遷移の計算結果が示されている。横軸はステップ数を表し、縦軸はＴであり、正の実数を表す。

図２４に示されているように、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合、各レプリカには各温度がほぼ等確率で設定され、各レプリカは温度空間をランダムウォークしていることが分かる。

図２５は、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合の各レプリカにおけるエネルギーの計算結果の一例を示す図である。図２５には、レプリカ番号＝１～８までの８つのレプリカのそれぞれにおける、エネルギー（Ｅ）の計算結果が示されている。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。

図２５に示されているように、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合、各レプリカのエネルギーは、特定の極小値に陥っても容易に脱出できていることが分かる。
図２６は、ボルツマン分布とＰｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合のレプリカ交換法の計算結果の一例を示す図である。また、図２７は、図２６の計算結果において、エネルギーが小さい領域を拡大した図である。図２６、図２７において、横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。

図２６から分かるように、ボルツマン分布を用いた場合（左）、エネルギーは特定の極小値に捕捉されるが、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合（右）、特定の極小値に陥っても容易に脱出できていることが分かる。また、サンプリングできるエネルギー領域も、ボルツマン分布を用いた場合よりも、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合のほうが広い。

さらに、図２７から分かるように、ボルツマン分布を用いた場合よりも、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合のほうが得られるエネルギーの最小値も小さい。したがって、最小値探索性能も、ボルツマン分布を用いた場合よりも、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合のほうが優れている。Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合、エネルギー差が大きい場合にも比較的大きな遷移確率をもつため、ボルツマン分布を用いる場合と異なり、温度の設定を粗くしても交換確率が稼げるからである。

図２８は、λアニール法の計算結果の一例を示す図である。図２８の左図は、λアニール法単独の場合の計算結果を示し、右図は、λアニール法とレプリカ交換法を組合せた場合の計算結果を示す。横軸はλを表し、縦軸はエネルギー（Ｅ）（＝Ｈ（λ））を表す。

使用した巡回セールスマン問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、計算条件は、式（２３）のｋ_１，ｋ_２＝１５０、Ｎ_Ｅ＝２００，０００（レプリカ交換法の場合はレプリカごと）、Ｎ_Ｉ＝１０、式（１４）のｍ_１＝は１．０２とした。温度は、λアニール法単独の場合、Ｔ＝１．０、レプリカ交換法と組合せた場合、８つのレプリカに対して、Ｔ＝１．０，３．０，５．０，１０．０，１５．０，２０．０，２５．０，３０．０が設定されるものとした。

λについてのアニールは、イジング型の量子コンピュータであれば、λを十分ゆっくり少しずつ増やしていけば基底状態が得られることが保証されている。しかし、量子論を扱わず、古典的なハミルトニアンの範囲内でエネルギーの最小値を求める最適化装置２０では、量子効果がないためその保証はない。図２８の左図の例では、λが大きくなるとサンプリング効率が悪化している。しかし、図２８の右図のように、λアニール法とレプリカ交換法を組合せることで、サンプリング効率を向上させることができる。これにより効率的にヒューリスティックに最小値の求解を行うことができる。

図２９は、デジタル回路を用いて並列試行を行う装置と、第２の実施の形態の最適化装置による計算結果の比較例を示す図である。横軸はエネルギーの小さい順序を示すランキングを表し、縦軸はエネルギー（Ｅ）を表す。

計算対象の問題は、上記と同様にｕｌｙｓｅｅｓ１６であり、式（２３）のｋ_１，ｋ_２は１５０である。
計算結果４０は、たとえば、非特許文献１に示されているようなデジタル回路を用いて並列試行を行う装置を用いた場合の計算結果である（ＳＡ法とダイナミックオフセット法を適用している）。試行数（前述のＮ_Ｅに対応）は、２億回である。

計算結果４１は、第２の実施の形態の最適化装置２０を用いた場合の計算結果である。なお、計算結果４１は、図５に示した処理フローに記載されているアニーリング処理やレプリカ交換処理が行われないときの計算結果を示している。Ｎ_Ｅ＝２００，０００、Ｎ_Ｉ＝１０である。最適化装置２０において、使用したＣＰＵ２１のクロック数は２．６０ＧＨｚであり、ＲＡＭ２２の容量は１２２８８ＭＢである。なお、プログラムはＣ＋＋で実装され、ＯｐｅｎＭＰやＭＰＩ（Message Passing Interface）による並列化は行われておらず、１ＣＰＵコアでの計算が行われた。

図２９に示されているように、デジタル回路を用いて並列試行を行う装置を用いた場合よりも、最適化装置２０を用いた場合のほうが、より小さい最小値が得られ、計算精度が優れていた。また、計算時間は、前者の場合、１３．５秒、後者の場合、１．６秒であり、後者のほうが計算速度についても優れていた。

次に、エネルギーが連続関数で表される場合の計算例を示す。
使用したエネルギー関数は、非特許文献３に記載されている以下の式（２４）で表される。

このエネルギー関数は、Ｅ_０≒５７．３２７６のときに、最小値０をとる関数である。式（２４）において、ｘ_ｉ（状態変数）の定義域は任意の実数である。このような実数で定義されたエネルギー関数についても、式（１４）で表されるＰｏｗｅｒ ｌａｗ型の遷移確率分布を用い、式（９）を用いたＭＣＭＣ計算を行うことで、効率的にエネルギー関数の最小値の探索が可能である。

図３０は、連続関数であるエネルギー関数に対する探索性能の比較例を示す図である。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）（式（２４）のエネルギー関数の値）を表す。図３０では、遷移確率分布としてボルツマン分布を用いた場合のサンプリング結果（左）、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合のサンプリング結果（右）が示されている。

計算条件は、Ｔ＝１．０、Ｎ_Ｅ＝１００，０００、Ｎ_Ｉ＝１０、式（１４）のｍ_１＝は１．０２である。
図３０から明らかなように、ボルツマン分布を遷移確率分布とした場合、一度、解が極小値に陥った場合、解が極小値から抜け出すことはない。一方、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合、極小値から抜け出すことができることが図３０から明らかである。また、Ｐｏｗｅｒ ｌａｗ型の遷移確率分布を用いた場合、図３０から明らかなように、一定回数同じ状態に留まるような状態を、少ない計算回数で多数求めることができる。すなわち、多数の極小値が求まるため、最小値以外にも極小値が重要な意味を持つ系では、効率的な解候補を求めることができる。そのような系の代表例として、有機化合物や生体分子などがある。

（変形例）
上記の説明では、図４に示した処理部３１は、ＣＰＵ２１が実行するプログラムモジュールを用いて実装できるものとして説明したが、処理部３１の各部を、ＡＳＩＣやＦＰＧＡなどの特定用途の電子回路により実現するようにしてもよい。

たとえば、ハミルトニアン計算部３１ｄは、積和演算回路を用いて式（１）を計算してもよい。ＭＣＭＣ計算実行部３１ｈについても、一様乱数と式（１４）～式（１６）に示されているようなｆ（ΔＥ）とを比較する比較回路などを用いて実現できる。さらに、エネルギー差計算部３１ｅを、非特許文献１などのように、各状態変数の値が変化することによるエネルギー差を並列に計算する電子回路により実現するようにしてもよい。最適化装置２０は、前述のようにボルツマン分布とは異なる遷移確率分布を用いることで、高効率で解を極小値から脱出させることができるが、上記のような電子回路をアクセラレータとして用いることで、解が深い極小値に陥った場合の脱出を促進できる。同様の理由により、ＧＰＧＰＵなどをアクセラレータとして処理部３１の各部の処理を実行するために用いてもよい。

なお、ＯｐｅｎＭＰやＭＰＩを用いた並列化や、ＡＶＸ（Advanced Vector eXtensions）２５６やＡＶＸ５１２など特定のＣＰＵ命令セットに対して特化したチューニングを行うことでより性能を向上できることが期待される。

なお、上記では、計算対象の問題として巡回セールスマン問題を例にして説明したが、金融工学など、他の分野の問題にも適用可能である。また、ジョブスケジューリング問題もイジングモデルを用いて定式化可能であるため、各種のスケジュール最適化問題（たとえば、病院内での関係者のスケジュールを最適化する問題など）に適用することもできる。さらに、ディープラーニング分野においても、制限ボルツマンマシンはイジングモデルを用いて定式化可能であるため、制限ボルツマンマシンの最適化にも使用できる。このため、人工知能分野にも適用可能である。

また、イジング型の評価関数の最小値の探索以外にも、状態変数が連続変数である評価関数の最小値の探索に対しても適用できるため、計算対象の問題も上記の分野に限定されるものではない。

また、上記では評価関数の最小値を探索するものとしたが、評価関数の符号などを反転することで、評価関数の最大値を探索する手法にも拡張できる。
なお、前述のように、上記の処理内容は、最適化装置２０にプログラムを実行させることで実現できる。

プログラムは、コンピュータ読み取り可能な記録媒体（たとえば、記録媒体２６ａ）に記録しておくことができる。記録媒体として、たとえば、磁気ディスク、光ディスク、光磁気ディスク、半導体メモリなどを使用できる。磁気ディスクには、ＦＤ及びＨＤＤが含まれる。光ディスクには、ＣＤ、ＣＤ－Ｒ（Recordable）／ＲＷ（Rewritable）、ＤＶＤ及びＤＶＤ－Ｒ／ＲＷが含まれる。プログラムは、可搬型の記録媒体に記録されて配布されることがある。その場合、可搬型の記録媒体から他の記録媒体（たとえば、ＨＤＤ２３）にプログラムをコピーして実行してもよい。

以上、実施の形態に基づき、本発明の最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムの一観点について説明してきたが、これらは一例にすぎず、上記の記載に限定されるものではない。

１０ 最適化装置
１１ 記憶部
１２ 処理部

Claims (12)

1. 問題を変換した評価関数に含まれる状態変数の値を記憶する記憶部と、
   現在の前記状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、前記状態変数の値が変化することによる前記評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により前記状態変数の値を更新する処理を繰り返すことで、前記評価関数の最小値の探索を行う処理部と、
   を有する最適化装置。
2. 前記遷移確率分布は、前記評価関数の値の変化と逆温度との積に１を加えた値のｍ（ｍ＞１）乗の逆数で表される、請求項１に記載の最適化装置。
3. 前記遷移確率分布は、前記評価関数の値の変化と逆温度との積の２乗に１を加えた値のｍ／２（ｍ＞１）乗の逆数で表される、請求項１に記載の最適化装置。
4. 前記遷移確率分布は、前記評価関数の値の変化と逆温度との積のｍ（ｍ＞０、ただしｍ≠１）乗にマイナスを掛けた値の指数関数で表される、請求項１に記載の最適化装置。
5. 前記処理部は、前記状態変数の値を更新するたびに前記評価関数の値を計算し、計算した前記評価関数の値のうち、小さい順に複数個の値を前記記憶部に記憶させる、請求項１乃至４の何れか一項に記載の最適化装置。
6. 前記処理部は、所定のサンプリング頻度で前記評価関数の値を出力する、請求項１乃至５の何れか一項に記載の最適化装置。
7. 前記処理部は、前記状態変数の値を更新する処理の回数の増加に伴って所定のスケジュールにより０から１まで増加するアニーリング変数を用い、１から前記アニーリング変数の値を引いた値と既知の第１のハミルトニアンとの積と、前記アニーリング変数と前記評価関数の値との積との和で表される第２のハミルトニアンを、前記アニーリング変数の各値について計算して計算結果を出力する、請求項１乃至６の何れか一項に記載の最適化装置。
8. 前記処理部は、前記アニーリング変数を０から１まで増加させる複数のスケジュールの中から、前記所定のスケジュールを選択する、請求項７に記載の最適化装置。
9. 前記処理部は、複数のレプリカのそれぞれに互いに異なる温度パラメータの値を設定し、設定された前記温度パラメータの値と前記遷移確率分布に基づいた交換確率にしたがって前記複数のレプリカ間で前記温度パラメータの値を交換するレプリカ交換法により前記最小値の探索を行う、請求項１乃至８の何れか一項に記載の最適化装置。
10. 前記処理部は、前記遷移確率分布に含まれるパラメータの値を変え、各パラメータの値における前記評価関数の値についての評価値を、前記評価関数の前記最小値の探索を行う際の計算時間よりも短い計算時間で計算し、計算結果を出力する、請求項１乃至９の何れか一項に記載の最適化装置。
11. 最適化装置の処理部が、
    記憶部に記憶されている問題を変換した評価関数に含まれる状態変数の値を取得し、
    現在の前記状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、前記状態変数の値が変化することによる前記評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により前記状態変数の値を更新する処理を繰り返すことで、前記評価関数の最小値の探索を行う、
    最適化装置の制御方法。
12. 記憶部に記憶されている問題を変換した評価関数に含まれる状態変数の値を取得し、
    現在の前記状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和を１に規格化できる関数で表されるとともに、前記状態変数の値が変化することによる前記評価関数の値の変化が正に大きいほど、ボルツマン分布よりも遷移確率が大きくなる遷移確率分布に基づいて、マルコフ連鎖モンテカルロ法により前記状態変数の値を更新する処理を繰り返すことで、前記評価関数の最小値の探索を行う、
    処理をコンピュータに実行させる最適化装置の制御プログラム。

Images