JP2022094510A - 最適化プログラム、最適化方法及び情報処理装置 - Google Patents

Abstract

【課題】レプリカ交換法における複数の温度パラメータの値の決定を容易にする。【解決手段】処理部１２は、評価関数に含まれる複数の状態変数の何れかの値が変化することによる評価関数の値の変化分と複数の温度パラメータの何れかの値とに基づいて得られる第１の遷移確率であって、温度パラメータの値の変化に対する評価関数の値の変化がボルツマン分布に基づく第２の遷移確率を用いた場合よりも緩やかになる第１の遷移確率にしたがって、複数の状態変数の何れかの値の更新を繰り返す更新処理を、複数のレプリカのそれぞれについて互いに独立に行うとともに、第１の遷移確率によって得られる確率分布の不変分布条件を満たす交換確率にしたがって、複数のレプリカの間で、複数のレプリカのそれぞれに設定された複数の温度パラメータの何れかの値、または複数のレプリカのそれぞれにおける複数の状態変数の値を交換する交換処理を繰り返す。【選択図】図１

Description

本発明は、最適化プログラム、最適化方法及び情報処理装置に関する。

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

イジングモデルを用いた最小値求解問題の計算手法として、マルコフ連鎖モンテカルロ法（以下ＭＣＭＣ法という）に基づき、疑似的な温度を示す温度パラメータを導入し高温から徐々に温度を下げる手法がある（たとえば非特許文献１，２，３参照）。この手法は、シミュレーテッド・アニーリング法（以下ＳＡ法と略す）と呼ばれる。ＳＡ法は最適解に到達することが理論的に保証されている方法である。しかし、ＳＡ法は、温度を対数の逆数にしたがって下げる手法であり、実用的ではない。そのため、実用上はそれよりも緩めたべき乗型の徐冷スケジュールが使われることが多いが、その場合、最小値への到達は保証されない。また、対数の逆数よりも高速な徐冷スケジュールを用いた場合、一度極小値にはまると極小値から抜け出せない欠点がある。

この欠点を考慮したアルゴリズムにレプリカ交換法がある（たとえば、特許文献１，２、非特許文献４参照）。レプリカ交換法では温度だけ異なる同じシミュレーションボックス（レプリカと呼ばれる）を多数用意し、一定頻度で交換条件を満たしたレプリカ間で温度を交換する。交換の結果、各レプリカは温度空間をランダムウォークする。その結果、高温領域では深いエネルギー極小値を抜け出すことができ、求解効率が上がる。

一方、レプリカ交換法にも欠点が存在する。レプリカ交換法は元々、物性物理学や計算化学分野で開発された手法であるため、レプリカ内部の遷移確率はボルツマン分布に基づいて指定される。そして、各レプリカの確率分布は、レプリカ交換を行ってもボルツマン分布に保たれるようにする。この条件は不変分布条件と呼ばれる。不変分布条件を課すのは統計物理学においては物理量の計算を前提とするためである。レプリカ交換法における問題点は系の自由度を増やしていくと、用いられるレプリカ数が自由度Ｎの関数としてＮの平方根に比例して増大する点である。

また、温度変化に対して評価関数の値（エネルギー）が急激に変化する相転移がある場合、レプリカ交換は機能しづらくなり、温度パラメータの刻み幅もより精密に設定しなければならず、計算者にとって負担になる。自由度が大きくなるほどこの傾向は強くなる。

この欠点を克服するアルゴリズムとしてマルチカノニカル法が提案されている（たとえば、非特許文献５参照）。マルチカノニカル法の場合、温度空間ではなく、エネルギー空間を等確率に訪問するようにアルゴリズムを組むことにより上記欠点を解決しようとするものであり、フラットヒストグラム法とも呼ばれる。しかし、マルチカノニカル法にも欠点がある。それはフラットヒストグラムを作成するために多くの予備計算を要する点である。さらにフラットヒストグラムが常に得られるとは限らず、計算量を多くしてもフラットヒストグラムにならない場合もある。そのため、予備計算をするために多くの労力を費やすことになる。

特開２０２０－０８６８２１号公報 特開２０１９－０７１１１９号公報 特開２０２０－０６４５３６号公報 特開２０１９－１９７３５５号公報 特開２０１８－０６７２００号公報

Sanroku Tsukamoto, Motomu Takatsu, Satoshi Matsubara and HirotakaTamura, "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, 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 Bernd A. Berg and Tarik Celik, "New Approach to Spin-Glass Simulations", Phys. Rev. Lett, Vol.69, No.15, October, 1992, p.2292 T. J. P. Penna, "Traveling salesman problem and Tsallisstatistics", Phys. Rev. E, Vol.51, No. 1, R1, January, 1995, p.51

上記のように、最適解の探索手法としてレプリカ交換法を用いた場合、サンプリング空間を広げ、より広範な解を求めつつ、求解効率を上げるための適切な温度パラメータの決定が難しい。

１つの側面では、レプリカ交換法における複数の温度パラメータの値の決定が容易な最適化プログラム、最適化方法及び情報処理装置を提供することを目的とする。

１つの実施態様では、最適化プログラムが提供される。最適化プログラムは、問題を変換した評価関数の情報を取得し、レプリカ交換法による最適解の求解処理に用いる互いに異なる複数の温度パラメータの値を決定し、前記複数の温度パラメータの値をそれぞれ複数のレプリカの何れかに１つずつ設定し、前記評価関数に含まれる複数の状態変数の何れかの値が変化することによる前記評価関数の値の変化分と前記複数の温度パラメータの何れかの値とに基づいて得られる第１の遷移確率であって、温度パラメータの値の変化に対する前記評価関数の値の変化がボルツマン分布に基づく第２の遷移確率を用いた場合よりも緩やかになる前記第１の遷移確率にしたがって、前記複数の状態変数の何れかの値の更新を繰り返す更新処理を、前記複数のレプリカのそれぞれについて互いに独立に行うとともに、前記第１の遷移確率によって得られる確率分布の不変分布条件を満たす交換確率にしたがって、前記複数のレプリカの間で、前記複数のレプリカのそれぞれに設定された前記複数の温度パラメータの何れかの値、または前記複数のレプリカのそれぞれにおける前記複数の状態変数の値を交換する交換処理を繰り返すことで、前記求解処理を実行する、処理をコンピュータに実行させる。

また、１つの実施態様では、最適化方法が提供される。
また、１つの実施態様では、情報処理装置が提供される。

１つの側面では、レプリカ交換法における複数の温度パラメータの値の決定が容易になる。

第１の実施の形態の情報処理装置の一例を示す図である。 ボルツマン分布に基づく遷移確率を用いた場合のエネルギーの温度変化の例を示す図である。 物理量Ａについての状態数とエネルギー区間分割との関係の例を示す図である。 ボルツマン分布に基づく遷移確率を用いた場合と、べき乗型の遷移確率を用いた場合のエネルギーの温度依存性の違いの例を示す図である。 べき乗型の遷移確率の指数とエネルギー（期待値）の温度依存性との関係の例を示す図である。 温度パラメータの各値が設定されるレプリカにおいて得られるエネルギー空間上での確率密度分布を示す図である。 レプリカ交換を行った場合に得られるエネルギー空間上の確率密度分布を示す図である。 レプリカ交換法による最大カット問題の計算例を示す図である。 レプリカ交換法における詳細つり合いの例を示す模式図である。 隣接温度パラメータ間のみによるレプリカ交換の例を示す図である。 エネルギー空間上の確率密度分布の頂点を与えるエネルギーと確率密度分布の標準偏差の例を示す図である。 低温領域での確率密度分布の振る舞いの例を示す図である。 第２の実施の形態の情報処理装置のハードウェアの一例を示す図である。 第２の実施の形態の情報処理装置の機能例を示すブロック図である。 第２の実施の形態の情報処理装置の一例の処理の流れを示すフローチャートである。 情報読込処理の一例の処理の流れを示すフローチャートである。 スピン初期化処理の一例の処理の流れを示すフローチャートである。 温度パラメータ計算処理の一例の処理の流れを示すフローチャートである。 確率密度の計算処理の一例の処理の流れを示すフローチャートである。 確率密度の更新処理の一例の流れを示すフローチャートである。 Ｅ_ｍｉｎとＥ_ｍａｘの更新処理の一例の流れを示すフローチャートである。 レプリカ交換処理の一例の流れを示すフローチャートである。 温度パラメータの値の交換の様子を示す図である。 レプリカに設定される温度パラメータの値の変化を示す図である。 ボルツマン分布に基づく遷移確率を用いた場合とべき乗型の遷移確率を用いた場合のレプリカ交換処理時のトンネル時間の比較結果の例を示す図である。

以下、発明を実施するための形態を、図面を参照しつつ説明する。
（第１の実施の形態）
図１は、第１の実施の形態の情報処理装置の一例を示す図である。

情報処理装置１０は、記憶部１１と処理部１２を有する。
記憶部１１は、問題を変換した評価関数（以下エネルギー関数という）の情報（以下問題情報という）を記憶する。また、記憶部１１は、評価関数に含まれる状態変数の値と状態変数の値に対応した評価関数の値（以下エネルギーという）の現在の最小値（最小エネルギー（Ｅ_Ｍｉｎ））などを記憶する。後述のように処理部１２は、レプリカ交換法によって問題の最適解（たとえば、エネルギー関数の最小値）を探索するものであるため、状態変数の値やエネルギーは、レプリカごとに記憶される。記憶部１１は、ＲＡＭ（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）などの特定用途の電子回路を含んでもよい。プロセッサは、ＲＡＭなどのメモリに記憶されたプログラムを実行する。たとえば、最適化プログラムが実行される。なお、複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」ということがある。

処理部１２は、レプリカ交換法により、たとえば、問題を変換したイジング型のエネルギー関数の最小値（または最小値が得られる状態変数の値の組合せ）を探索する。
イジング型のエネルギー関数（Ｈ（｛ｘ｝））は、たとえば、以下の式（１）で定義される。

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

重み係数（Ｗ_ｉｊ）、バイアス係数（ｂ_ｉ）、定数（Ｃ）は、問題情報として、記憶部１１に記憶される。
ここで、以下のように、記号を定義する。まず、Ｎ個の状態変数の組合せによって得られる状態は、式（１）から離散有限個存在する。状態の総数をＭとすると、Ｍ＝２^Ｎである。また、エネルギーの値が異なる状態数をＭ_Ｅとすると、Ｍ_Ｅ≦２^Ｎである。エネルギーの値を小さい順にＥ_０，Ｅ_１，…，Ｅ_ｋ，…，Ｅ_ＭＥと記述する。また、Ｅ_ｋを与えるＮ個のｘ_ｉの組を｛ｘ_ｋ｝と記述する。したがって、Ｈ（｛ｘ_ｋ｝）＝Ｅ_ｋとなる。なお、このとき、縮重度をＭ_ｋとする。同じエネルギーを与える状態変数の組合せは複数存在することが多い。

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

上記のようなエネルギー関数の最小値を探索する際、処理部１２は、レプリカ交換法を用いる。
レプリカ交換法では、式（１）で定義されるエネルギー関数の値をそれぞれ計算する複数のレプリカが用意される。レプリカ番号＝ｉのレプリカに設定される温度パラメータをＴ_ｉとする。各レプリカでは固定の温度パラメータ（前述の疑似的な温度）の値を用いて一定回数のＭＣＭＣ計算が行われる。そして、その一定回数ごとに、所定の交換確率に基づいて、レプリカ間における温度パラメータの値の交換が行われる。なお、温度パラメータの値の交換の代わりに、状態（Ｎ個の状態変数の値）を交換してもよい。

ＭＣＭＣ計算において、処理部１２は、所定の遷移確率にしたがって状態遷移を発生させる。このとき、エネルギーが大きくなる状態遷移についても一定確率で許容される。これはメトロポリス法として知られている。従来のようにボルツマン分布を用いた場合、メトロポリス法における遷移確率は、Ｐ_ｉ→ｊ＝ｍｉｎ（１，ｅｘｐ（－βΔＥ））、と表せる。Ｐ_ｉ→ｊは、状態遷移前の状態ｉと、状態ｉにおける状態変数からランダムに１つの状態変数を選び反転させることによる状態遷移後の状態ｊの間の遷移確率である。βは、各レプリカに設定される温度パラメータの値（ここでは温度の逆数（逆温度と呼ばれる））であり、レプリカ番号＝ｉのレプリカにおける逆温度は、β_ｉ＝１／Ｔ_ｉである。

レプリカ交換法における交換確率は、従来のようにボルツマン分布に基づく遷移確率を用いた場合、以下の式（３）で表せる。

式（３）において、Ｐ_Ａ（ｔ）は、レプリカ交換前の状態Ａが実現する確率である。Ｐ_Ｂ（ｔ）は、レプリカ交換後の状態Ｂが実現する確率である。Δβは、レプリカ番号＝ｊのレプリカに設定されている逆温度（β_ｊ）とレプリカ番号＝ｉのレプリカに設定されている逆温度（β_ｉ）との差であり、Δβ＝β_ｊ－β_ｉである。ΔＥは、レプリカ番号＝ｊのレプリカのエネルギー（Ｅ_ｊ）とレプリカ番号＝ｉのレプリカのエネルギー（Ｅ_ｉ）との差であり、ΔＥ＝Ｅ_ｊ－Ｅ_ｉである。

つまり、レプリカ交換法では、レプリカ同士の逆温度差とエネルギー差の積で決まる確率で温度パラメータの値が交換される。温度パラメータの値を交換することにより、レプリカ交換法ではそれぞれのレプリカが温度空間をランダムウォークするようになる。その結果、高温領域を経由することができるようになり、状態が、エネルギーが深い極小値をもつ極小状態にはまったとしても容易にそこから脱出できる。

しかし、レプリカ交換法において、ボルツマン分布に基づいて交換確率を式（３）とした場合、２つのレプリカ同士のエネルギーが離れるほど交換確率は指数関数的に下がってしまう。

このため、温度パラメータの設定には制限があり、温度パラメータの刻み幅は概算でＮの平方根分の１程度、レプリカ数は概算でＮの平方根程度必要になる。これは自由度を増やしたときにレプリカ数が増大することを意味する。

また、温度パラメータの設定について、以下のような問題がある。
図２は、ボルツマン分布に基づく遷移確率を用いた場合のエネルギーの温度変化の例を示す図である。横軸は温度パラメータ（Ｔ）、縦軸はエネルギー（Ｅ）を表している。

図２では、ある問題を計算した場合に得られるエネルギーの温度変化の例が示されている。この例の場合、１０＜Ｔ＜１００の範囲において、エネルギーが温度パラメータの値に対して急激に増大していることが分かる。これは相転移にまつわる現象である。有限自由度の系では相転移は起きないことが知られているが、系の自由度が増えるほど、エネルギーの温度変化が急激に起こるようになり、相転移の特徴が際立っていく。

このような相転移が起こると、レプリカ交換が効率的に機能しなくなる。つまり、相転移が生じる温度パラメータの値（以下転移温度という）以下の領域でのみレプリカ交換が起こるグループと、転移温度より大きい領域でのみレプリカ交換が起こるグループに分かれてしまう。このため、自由度が大きくなるほどレプリカ交換を機能させるためには、転移温度近傍でのレプリカの温度刻みをより慎重に選択しなければならなくなる。

この原因になっているのは式（３）のような交換確率を用いるためである。式（３）は、レプリカ交換を行っても各レプリカがしたがう確率分布がボルツマン分布のまま保たれる条件で定式化したものである。このような条件は、不変分布条件と呼ばれている。具体的には状態Ａをレプリカ番号＝ｉのレプリカが（｛ｘ_ｉ｝，β_ｉ）かつ、レプリカ番号＝ｊのレプリカが（｛ｘ_ｊ｝，β_ｊ）である状態とする。そして、状態Ｂをレプリカ番号＝ｉのレプリカが（｛ｘ_ｊ｝，β_ｉ）かつ、レプリカ番号＝ｊのレプリカが（｛ｘ_ｉ｝，β_ｊ）である状態とする。交換の前後で両者の確率分布が変わらないとして任意の２つのレプリカ間で定義されるレプリカ交換が定常状態に達する条件は、以下の式（４）で与えられる。

式（４）において、π（｛ｘ_ｊ｝，β_ｊ）は系がしたがう統計分布である。π（｛ｘ_ｊ｝，β_ｊ）としてボルツマン分布を用いると、式（３）の交換確率の式が得られる。
磁性体などの物理系についての問題を扱う場合、ボルツマン分布を使わざるを得ない。なぜならば、熱平衡状態がｅｘｐ（－βＥ）の統計分布にしたがうからである。また、ヘルムホルツ自由エネルギーやギブス自由エネルギーなど、系の熱力学的挙動にエントロピー効果が重要な役割を果たす場合は、エントロピーの効果を適切に取り入れるため、使用可能な確率分布にも制約がかかる。

しかしながら、イジングモデルで表現される式（１）のようなエネルギー関数の最小値を求める問題は、単に関数の最小値求解問題とみなせばよいため、ボルツマン分布に縛られなくてもよい。ボルツマン分布により相転移が引き起こされ、相転移がレプリカ交換の効率を阻害するのであれば、相転移が起きない確率分布を用いればよい。

さらにマルチカノニカル法のように予備計算が面倒であるならば、予備計算が楽になるような計算法を用いればよい。マルコフ連鎖が既約（irreducible）でありさえすれば、系の最小値にたどり着くことが確率論的に保証される。

そこで、処理部１２は、ボルツマン分布使用時に見られる相転移現象を回避するため、ボルツマン分布ではない確率分布を用いる。
以下では、遷移確率を上記のようにＰ_ｉ→ｊ＝ｍｉｎ（１，ｅｘｐ（－βΔＥ））と定義するのではなく、以下の式（５）のように定義する。

式（５）において、遷移確率ｆ（ΔＥ_ｉｊ）は任意の関数であるが、有限確定し、ｆ（ΔＥ_ｉｊ）＜∞を満たすものとする。また、境界条件として、以下の式（６）が要請される。

さらに、現在の状態変数の値により表される現在の状態から複数の異なる状態のそれぞれに遷移する確率の和が有限確定するものとし、以下の式（７）の条件が満たされるものとする。

なお、式（７）においてＢ_ｉは、たとえば、Ｂ_ｉ＝１である。ただし、Ｂ_ｉ＝１に限定されるわけではない。
次に、式（５）で定義される遷移確率の下での不変分布条件を満たすレプリカ交換を定義する。式（５）で定義される遷移確率の下で定まる確率分布π（｛ｘ_ｉ｝，ｔ）は、以下の式（８）のマスター方程式によって決定される。

π（｛ｘ_ｉ｝，ｔ）は、式（５）の遷移確率の定義より一意に存在する。なぜならば、全てのΔＥ_ｉｊに対する遷移確率が０ではなく、状態を変えない遷移であるΔＥ_ｉｉ＝０に対する遷移確率も０ではないからである。式（８）において定常状態は、以下の式（９）のように与えられる。

このため、定常状態が満たすべき方程式は、以下の式（１０）で表せる。

定常状態では、確率分布は時刻ｔに依存しないためｔを省いてπ（｛ｘ_ｉ｝）と表記できる。式（５）の遷移確率により導入される確率分布では、詳細つり合いの原理が満たされない場合がある。詳細つり合いの原理とは、式（１０）における総和記号の個々の項が０になるというものである。詳細つり合いが成り立つ場合、Ｐ_ｍ→ｊπ（｛ｘ_ｍ｝）－Ｐ_ｊ→ｍπ（｛ｘ_ｊ｝）＝０が成立するため、以下の式（１１）が成立する。

式（１１）から、確率分布をボルツマン分布型としてｅｘｐ（－βΔＥ）を用いると、遷移確率は、前述のＰ_ｉ→ｊ＝ｍｉｎ（１，ｅｘｐ（－βΔＥ））となる。式（５）の遷移確率の下では詳細つり合いの原理は成立しないが、異なる２つの状態変数の組｛ｘ_ｋ｝と｛ｘ_ｌ｝（ｋ≠ｌ）に対して、同一のエネルギー（Ｅ_０）が与えられるとき、確率分布関数は同じ値になる。つまり、π（｛ｘ_ｋ｝）＝π（｛ｘ_ｌ｝）が成立する。

これは以下の理由による。２つの異なる状態変数の組｛ｘ_ｋ｝と｛ｘ_ｌ｝（ｋ≠ｌ）に対してＥ_０となるため、Ｅ_０＝Ｈ（｛ｘ_ｋ｝）＝Ｈ（｛ｘ_ｌ｝）となる。ここで、式（１０）から、以下の式（１２）が得られる。

式（１２）において、Ｐ_ｌ→ｋ＝１、Ｐ_ｋ→ｋ＝１であるのは明らかである。同様に、以下の式（１３）が得られる。

遷移先の｛ｘ_ｋ｝と｛ｘ_ｌ｝はエネルギーの値が同じである。したがって、遷移確率も同じであり、以下の式（１４）が成立する。

式（１３）から式（１２）を差し引くと、π（｛ｘ_ｌ｝）Ｂ_ｌ－π（｛ｘ_ｋ｝）Ｂ_ｋ＝０が成立する。ここで、２つの異なる状態変数の組｛ｘ_ｋ｝と｛ｘ_ｌ｝（ｋ≠ｌ）に対してＥ_０となるため、Ｂ_ｋ＝Ｂ_ｌが成立する。そしてＢ_ｋ、Ｂ_ｌは０ではない。したがって、π（｛ｘ_ｌ｝）＝π（｛ｘ_ｋ｝）となる。

この条件は強力な要請であり、式（５）の遷移確率をエネルギー差だけで定義した場合、エネルギーの値が同じ状態は全て同じ確率で実現する。そして、この条件はエネルギーに限らず成立する。式（５）の遷移確率を任意の物理量Ａの差のみで定義したとき、物理量Ａが同じになる微視的状態は、詳細つり合い原理の成立の有無に関わらず全て同じ確率で実現する。ただし、これは確率分布が存在することは保証するが、解析的に記述できることは保証しない。また、べき乗型の遷移確率を用いたからといって、確率分布がべき乗型になるわけではない。

ボルツマン分布の場合、エネルギーの値のみで実現確率が定義され、エネルギーが同じ微視的状態の実現確率が同じになることを考慮すれば、この条件はボルツマン分布を特殊な場合として含む一般化を行っていることになる。なお、この条件は式（１）で表されるようなイジングモデルに限らず一般の離散有限状態系について成立する。

次にエネルギー空間での確率分布について説明する。
エネルギーの値がＥ_ｊになる状態変数の組｛ｘ_ｋ｝は複数存在する場合があるが、それらをまとめて｛ｘ_ｋ ^（ｊ）｝と表記する。このとき、状態の総数Ｍは以下の式（１５）で表せる。

図３は、物理量Ａについての状態数とエネルギー区間分割との関係の例を示す図である。横軸はＡを表し、縦軸は頻度を表す。
物理量Ａの各値（Ａ_１，Ａ_２，…，Ａ_ｊ，…Ａ_ＭＡ）に対して、同じ物理量Ａの状態変数の組が複数ある。

同じエネルギーに属する微視的状態（図３の状態変数の組｛ｘ_１｝、｛ｘ_２｝、｛ｘ_３｝など）への遷移はすべて同じ遷移確率になる。このことから、ｋ番目の微視的状態（｛ｘ_ｋ｝）がＥ_ｊに属するとして、Ｅ_ｊからＥ_ｍ’の状態への遷移確率をＰ_ｊ→ｍ’ ^（Ｅ）と表記することにすると、以下の式（１６）が得られる。

同様に、以下の式（１７）が得られる。

したがって、マスター方程式の定常解は以下の式（１８）のようになる。

ここで、ｎ^（Ｅ）（Ｅ_ｊ）＝Ｍ_ｊπ（Ｅ_ｊ）となる確率分布ｎ^（Ｅ）（Ｅ_ｊ）を導入すると、式（１８）は、以下の式（１９）のように表せる。

式（１９）の両辺にＭ_ｊを掛けると、以下の式（２０）が得られる。

ここで、遷移確率を定義し直し、以下の式（２１）のように定義する。

このとき、式（２０）は以下の式（２２）のように表せる。

したがって、エネルギーの値の縮退度の分だけ、エネルギー空間での遷移確率は増大する。実現確率もエネルギーの値の縮退の分だけ増大する。縮退度は状態数であるため、Ｎ＞＞１の場合、Ｗａｎｇ－Ｌａｎｄａｕ法などを用いれば近似的に求めることができ、微視的状態の実現確率も求めることができる。

次に、式（５）の遷移確率における確率分布を不変分布とするレプリカ交換について説明する。以下ではレプリカ数をＮ_{ｒｅｐｌｉｃａ}とし、全てのレプリカからなる系が状態｛｛ｘ_ｉ｝，β_ｉ，Ｅ_ｉ｝（ｉ＝１，２，…，Ｎ_{ｒｅｐｌｉｃａ}）で与えられる状態を状態Ａと呼び、状態Ａが実現される確率をＰ_Ａと記述することにする。このとき、Ｐ_Ａは、以下の式（２３）で表せる。

次に、レプリカ番号＝ｉのレプリカとレプリカ番号＝ｊのレプリカとの間で、レプリカ交換により温度パラメータの値を交換した状態を状態Ｂと呼び、状態Ｂが実現される確率をＰ_Ｂと記述することにする。このとき、Ｐ_Ｂは、以下の式（２４）で表せる。

状態Ａと状態Ｂの２状態がしたがうマスター方程式は、以下の式（２５）となる。

式（２５）において、Ｐ_Ａ→Ｂは状態Ａから状態Ｂへの遷移確率を示す。式（２５）のマスター方程式の定常状態は、Ｐ_Ａ（ｔ）Ｐ_Ａ→Ｂ－Ｐ_Ｂ（ｔ）Ｐ_Ｂ→Ａ＝０、という式により与えられる。

したがって、以下の式（２６）が得られる。

これにより、レプリカ交換の交換確率Ｐ_ｅｘは以下の式（２７）のように定義できる。

式（２６）の条件は、不変分布条件と呼ばれる。なぜならば、レプリカ交換が、確率分布の関数形を保存するような束縛条件が課されて行われるからである。
Ｐ_ｅｘを定義するための確率分布の解析形を求めることは困難であるが、この問題は以下のように解決できる。

前述のｎ^（Ｅ）（Ｅ_ｊ）＝Ｍ_ｊπ（Ｅ_ｊ）を用いて、式（２６）をエネルギー空間の表式に直すと、以下の式（２８）のようになる。

したがって、エネルギー空間表示をしたときの確率密度の比から、式（５）で表せる遷移確率によって得られる確率分布を不変分布にするようなレプリカ交換の交換確率を定義できる。これにより、各レプリカがしたがう確率分布がπ（｛ｘ_ｊ｝，β_ｊ，Ｅ_ｊ）にしたがうことを強制するレプリカ交換が実現される。つまり、この方式では、個々のレプリカ内部では詳細つり合いが要求されないが、レプリカ交換において詳細つり合いが要求される。式（２８）において、確率分布をボルツマン分布型にとると式（３）が得られる。

したがって、ボルツマン分布を用いたレプリカ交換ではメトロポリス法の式がたまたまよく似た式になるが一般には異なる。レプリカ内部での詳細つり合いの成立の有無と、レプリカ交換における詳細つり合いは理論的に別の起源をもつ。

そして、不変分布条件が課されるレプリカ交換は各レプリカが同一の確率分布を保持する条件であるため、この条件は最小値求解問題では一般には必要はない。もちろん、個々のレプリカ内の確率分布を、全体分布を構成するための材料とみなすのであれば、材料である確率分布は不変であったほうが制御はしやすい。しかし、一般には全てのレプリカを含む全系において既約でありさえすれば、解への到達可能性は保証される。

したがって、レプリカ内の遷移確率とレプリカ間の交換確率は別々に設定してよい。そのための条件は２つである。条件の１つ目はあるレプリカと別のレプリカとの交換確率が０ではないこと、条件の２つ目はあるレプリカから自分自身であるレプリカへの遷移が０ではないこと、である。この２つの条件が守られれば、計算者は計算に都合のよいレプリカ交換の交換確率を定義してよい。

図１の処理部１２は、以下の式（２９）の交換確率（Ｐ_ｅｘ）にしたがって、レプリカ交換を行う。

式（２９）のＰ_ｅｘを用いることで、不変分布条件が満たされる。
図１には、処理部１２による最適化方法の一例の処理の流れが示されている。
処理部１２は、記憶部１１から問題情報を取得する（ステップＳ１）。なお、処理部１２は、記憶部１１から式（５）のｆ（ΔＥ_ｉｊ）（遷移確率）の情報や、レプリカ交換による求解処理の終了条件となる計算回数など記憶部１１から取得してもよい。

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

その後、処理部１２は、レプリカ交換法による最適解の求解処理に用いる互いに異なる複数の温度パラメータの値を決定するとともに、複数の温度パラメータの値をそれぞれ複数のレプリカの何れかに１つずつ設定する（ステップＳ３）。

前述のように、ボルツマン分布に基づく遷移確率を用いた場合、図２のように、温度パラメータの値の変化に対するエネルギーの変化が急な部分があり、複数のレプリカのそれぞれに設定する温度パラメータを決定することが難しい。そのため、式（５）のｆ（ΔＥ_ｉｊ）として、温度パラメータの値の変化に対するエネルギー関数の値の変化がボルツマン分布よりも緩やかになる遷移確率が用いられる。これにより、温度パラメータの決定が容易になる。なお、温度パラメータの決定方法については後述の第２の実施の形態において説明する。

そして、処理部１２は、レプリカ交換法による求解処理を行う（ステップＳ４）。
処理部１２は、上記の遷移確率にしたがって、複数の状態変数の何れかの値の更新を繰り返す更新処理（一定回数のＭＣＭＣ計算）を、複数のレプリカのそれぞれについて互いに独立に行う。図１には、温度パラメータの値の変化に対するエネルギーの変化がボルツマン分布よりも緩やかになる遷移確率の例として、べき乗型の遷移確率（ｆ（ΔＥ_ｉｊ）＝１／（１＋βΔＥ）^ｍ（ｍ＞１）が示されている。

処理部１２は、上記一定回数ごとに、式（２９）のＰ_ｅｘにしたがって、複数のレプリカの間で、複数のレプリカのそれぞれに設定された複数の温度パラメータの何れかの値を交換する処理を繰り返す。なお、処理部１２は、温度パラメータの値の交換の代わりに、状態（Ｎ個の状態変数の値）を交換してもよい。

なお、式（２９）の確率密度（ｎ（β_ｉ，Ｅ_ｊ）など）は、温度パラメータの値ごとに独立なサンプリング計算を行うことで比較的容易に得られる。確率密度の計算方法の例やレプリカ交換法による求解処理のより詳細な例については、第２の実施の形態において説明する。

処理部１２は、たとえば、各レプリカにおいてＭＣＭＣ計算が行われるたびに、エネルギーを計算し、記憶部１１に記憶される現在のＥ_Ｍｉｎよりも低いエネルギーが得られた場合には、Ｅ_Ｍｉｎを更新する。そして、処理部１２は、所定の回数のレプリカ交換が終了した時点でのＥ_Ｍｉｎを計算結果として、たとえば、外部装置（外部のコンピュータ、記憶媒体、表示装置など）に出力し（ステップＳ５）、処理を終える。

なお、処理部１２は、Ｅ_Ｍｉｎが得られたときのｘ_１～ｘ_Ｎの値を記憶部１１に記憶して、最後に記憶したｘ_１～ｘ_Ｎの値をＥ_Ｍｉｎとともに出力してもよい。
以上のような第１の実施の形態の情報処理装置１０及び最適化方法によれば、各レプリカについて、温度パラメータの値の変化に対する評価関数の値の変化がボルツマン分布に基づく遷移確率よりも緩やかになる遷移確率を用いてＭＣＭＣ計算が行われる。これにより、ボルツマン分布に基づく遷移確率を用いた場合に生じる相転移が抑制されるため温度パラメータの決定が容易になる。

なお、ボルツマン分布に基づく遷移確率とは異なる遷移確率を用いた場合、不変分布条件は自明ではないが、上記のように、不変分布条件を満たす交換確率によりレプリカ交換を行うことで、レプリカ交換前後において確率分布が同じとなる。これにより、計算の安定性が増し、計算が制御しやすくなる。つまり、エネルギー空間を安定してサンプリング可能になり、その結果、求解効率が安定する。

（第２の実施の形態）
以下に示す第２の実施の形態では、温度パラメータの値の変化に対する評価関数の値の変化がボルツマン分布に基づく遷移確率よりも緩やかになる遷移確率として、以下の式（３０）で表されるべき乗型の遷移確率が用いられる。

図４は、ボルツマン分布に基づく遷移確率を用いた場合と、べき乗型の遷移確率を用いた場合のエネルギーの温度依存性の違いの例を示す図である。横軸は温度パラメータ（Ｔ）、縦軸はエネルギー（Ｅ）を表す。なお、用いたべき乗型の遷移確率は、式（３０）において、ｍ＝１．００１としたものである。

なお、エネルギーの期待値〈Ｅ〉は、以下の式（３１）で表せる。

Ｅ（Ｘ）は、式（１）のＨ（｛ｘ｝）であり、Ｐ（Ｘ）はある状態Ｘの確率分布であり、Ｎ_ｄａｔａは、サンプリングで取得したデータ点数である。
図４は、温度パラメータの値ごとに平衡状態に達してからサンプリングを開始し、十分大きい値であるＮ_ｄａｔａのＥ_ｉを取得することによって得られたものである。

図４のように、低温領域と高温の極限では、ボルツマン分布に基づく遷移確率を用いた場合と、べき乗型の遷移確率を用いた場合のエネルギーの温度依存性はほぼ等しい。しかし、中間の温度領域では顕著な違いが現れている。ボルツマン分布に基づく遷移確率を用いた場合、相転移に対応する温度パラメータの値の近傍でエネルギーが急激に変化している。これは自由度が増えるほど急激になり、λ点的になっていく。一方、べき乗型の遷移確率を用いた場合、エネルギーの増加が、ボルツマン分布に基づく遷移確率を用いた場合よりも緩やかになっている。

図５は、べき乗型の遷移確率の指数とエネルギー（期待値）の温度依存性との関係の例を示す図である。横軸は温度パラメータ（Ｔ）、縦軸はエネルギー（Ｅ）を表す。なお、比較のためにボルツマン分布に基づく遷移確率を用いた場合のエネルギーの温度依存性（“Ｂｏｌｚ”と表記されている）も示されている。

べき乗型の遷移確率の指数（ｍ）を増やしていくと徐々にエネルギーの温度に対する増加率が増大していく。つまり、ｍが増えるほど、より顕著に相転移らしき現象が現れてくる。したがって、たとえば、ｍ＞４の場合は、レプリカ交換時の温度パラメータの刻み幅をより慎重に選ぶことになる。温度パラメータの値の決定をより容易にするためには、たとえば、１＜ｍ≦４であることが望ましい。

図６は、温度パラメータの各値が設定されるレプリカにおいて得られるエネルギー空間上での確率密度分布を示す図である。横軸はエネルギー（Ｅ）、縦軸は確率密度ｎ（Ｅ）を表す。

なお、レプリカ交換は実施せず、各レプリカの温度パラメータの値は固定されている。また、図６の４つのケースの全ての計算において同じ温度列（同じ温度パラメータ群）が用いられている。なお、図６の計算例では、最低温度を示す温度パラメータ（Ｔ_ｍｉｎ）については、Ｔ_ｍｉｎ=１．０、最高温度を示す温度パラメータ（Ｔ_ｍａｘ）については、Ｔ_ｍａｘ＝１００、レプリカ数は２６としている。

図６では、４つのケースのエネルギー空間上の確率密度分布、すなわち、ボルツマン分布に基づく遷移確率を用いた場合に得られる確率密度分布と、ｍ＝１．００１，２，３としたべき乗型の遷移確率を用いた場合に得られる確率密度分布の例が示されている。

図６のように、相転移点近傍の中間のエネルギー状態（Ｅ＝－４０００～－１００００程度の領域）において、確率密度分布の頂点の間隔が、ボルツマン分布に基づく遷移確率を用いた場合、べき乗型の遷移確率分布を用いた場合よりも大きくなる傾向にある。これは温度パラメータの関数としてエネルギーが急激に変化することに対応する。

べき乗型の遷移確率を用いた場合は上記中間のエネルギー状態のエネルギーにおける頂点の間隔は、ボルツマン分布に基づく遷移確率を用いた場合よりも小さくなる傾向にある。このため、上記中間のエネルギー状態におけるサンプリングはボルツマン分布に基づく遷移確率を用いた場合よりも、サンプリングしやすい。

そして、べき乗型の遷移確率の指数であるｍが大きくなるにつれて、上記中間のエネルギー状態における確率密度分布の頂点間の間隔は大きくなる傾向にある。
つまり、隣接する温度パラメータの値が設定される２つのレプリカにおいて得られる確率密度分布の頂点間の間隔は、温度パラメータの値に対して鈍感になる。このため、ｍを小さくするほど温度パラメータの設定がより容易になる。

ただし、ボルツマン分布に基づく遷移確率とべき乗型の遷移確率が作り出す確率分布は本質的に異なる。そのため、その影響が低温領域に出ている。図６の計算例の場合、ボルツマン分布に基づく遷移確率を用いた場合にサンプリングがしやすいように温度パラメータを決定したため、べき乗型の遷移確率では低温領域でサンプリング能力が低下していることが分かる。これはべき乗型の遷移確率を用いた場合、ボルツマン分布に基づく遷移確率を用いた場合に十分低温と考えられる領域も十分低温でないからである。この問題は最低温度を示す温度パラメータ（Ｔ_ｍｉｎ）を、ボルツマン分布に基づく遷移確率を用いた場合よりも小さく取ることで解決できる。

図７は、レプリカ交換を行った場合に得られるエネルギー空間上の確率密度分布を示す図である。図７の計算例では、図６の計算時と同様の温度パラメータの設定条件が用いられている。図７の左側はサンプリングを行った全エネルギー領域について確率分布関数を数値的に計算したものである。図７の右側は図７の左側において、Ｅ＝－１２５００～－１２０００までを拡大したものである。横軸はエネルギー（Ｅ）、縦軸は確率密度（Ｐ（Ｅ））を表す。なお、用いたべき乗型の遷移確率は、式（３０）において、ｍ＝３としたものである。また、比較のためにボルツマン分布に基づく遷移確率を用いた場合の確率密度分布（“Ｂｏｌｚ”と表記されている）も示されている。

図７の右側のように、低エネルギー側ではボルツマン分布に基づく遷移確率を用いたほうが、べき乗型の遷移確率を用いた場合よりも確率密度が大きくなっている。一方で、図７の左側のように、高エネルギー側ではべき乗型の遷移確率を用いたほうが確率密度は大きくなっている。全体的にはべき乗型の遷移確率を用いた場合、エネルギー空間での確率密度の振る舞いは、ボルツマン分布に基づく遷移確率を用いた場合よりも平坦になっている。マルチカノニカル法の場合、数値的にエネルギー空間上で意図的に全エネルギー領域を等確率で訪問するようなフラットヒストグラムを作るが、マルチカノニカル法ほど平坦にはなっていない。

この特徴は基底状態を探索する場合に有利になる局面が存在することを意味する。それは効率的に記憶を忘却することが効率的な探索に繋がると考えられているからである。ここで述べた記憶の忘却とは、高温状態に遷移することによって特定のエネルギーランドスケープ上の構造の影響から脱するという意味である。レプリカ交換法の長所は高温側に状態を遷移させることで深いエネルギー極小構造からも効率的に抜け出せる点にある。しかし、系の自由度が大きくなると、相転移点近傍におけるエネルギー変化が急激になってしまい、効率的に高温領域に遷移させることが難しくなってしまう。逆に、効率的に高温領域に遷移させることができるということは、より広い解空間をサンプリングできることに繋がる。

これまでに説明した結果を使うと、べき乗型の遷移確率を用いた場合、レプリカ数を減らせることが分かる。
図８は、レプリカ交換法による最大カット問題の計算例を示す図である。採用された問題は、最適解が知られているＧ４３と呼ばれる問題である。横軸はレプリカ数、縦軸は求解確率（％）を表す。なお、用いたべき乗型の遷移確率は、式（３０）において、ｍ＝１．００１としたものである。また、比較のためにボルツマン分布に基づく遷移確率を用いた場合の計算例も示されている。

なお、図８は、各レプリカについて１０万回の計算が行われ、レプリカ交換頻度を１０回に一度としたものである。計算が終了したのち、基底状態の解（最適解）に到達したレプリカの数を数え、そこから求解確率が計算されている。たとえば、２６個のレプリカを用いて計算が行われた場合、２６個の全てのレプリカが最適解に到達したら求解確率は１００％である。

図７から分かるように、レプリカ数が少なくなると求解確率も減っていく。しかし、べき乗型の遷移確率を用いた場合、ボルツマン分布に基づく遷移確率を用いたときよりも少ないレプリカ数でも最適解が得られている。

これは先ほど述べたように、べき乗型の遷移確率を用いることで、エネルギーの温度依存性が温度パラメータの値に関してボルツマン分布に基づく遷移確率を用いた場合よりも鈍感になることによる効果である。つまり、温度パラメータの値を細かく取らなくてもよいため、レプリカ数もそれに伴い少なくすることができる。

以下この現象を、数式を使って説明する。いま、レプリカ交換が理想的に行われたとして、エネルギー空間を満遍なく平等に訪問するような確率分布であるρ（Ｅ）が得られたとする。ρ（Ｅ）は全てのレプリカによって得られる確率分布とする。このときρ（Ｅ）をボルツマン分布であるρ_Ｂの和として表現すると、以下の式（３２）のように表せる。

同様に、ρ（Ｅ）を、式（５）の遷移確率によって得られる確率分布であるρ_Ｐの和として表現すると、以下の式（３３）のように表せる。

つまり、効率的な計算をするためのρ（Ｅ）を、異なる２つの分布関数で構成する問題に帰着する。この場合、ボルツマン分布に基づく遷移確率を用いると、図５からも分かる通り、相転移点近傍でエネルギーが急激に変化する。したがって、式（３２）において、全てのエネルギー領域を等しくサンプリングしようとすれば、相転移点近傍において温度パラメータを多く取ることになる。そのため、式（３２）を構成する和の項数は多くなる。

一方、べき乗型の遷移確率においては図５からも分かる通り、相転移に対応する相転移点が明確に現れず、温度パラメータを多く取らずとも式（３３）から全体の望ましい確率分布を作ることができる。

次にレプリカ数を最小にするための条件について述べる。これは最小値求解には必須ではないが、レプリカ数を最小にするためには必須の条件である。その条件とは、レプリカ番号＝ｉのレプリカが作り出すエネルギー空間上の確率密度分布は、レプリカ番号＝ｉ－１及びレプリカ番号＝ｉ＋１のレプリカが作り出す確率密度分布とのみ重なり合いをもつという条件である。

図９は、レプリカ交換法における詳細つり合いの例を示す模式図である。
図９には、レプリカ番号＝１～４までの４つのレプリカについての詳細つり合いの様子が示されている。４つのレプリカのうち、任意の２つのレプリカの間でレプリカ交換が可能となっている。

レプリカ交換法においては、任意の２つの交換対象のレプリカに対して詳細つり合いの原理が課される。これは、式（２６）から要請される。詳細つり合いの原理は全ての任意の２つのレプリカについて要請される。そのため、一般には１回のレプリカ交換ではＮ_{ｒｅｐｌｉｃａ}（Ｎ_{ｒｅｐｌｉｃａ}－１）／２通りの組合せが生じる。しかし、レプリカ交換では温度パラメータの設定を適切に行うことでレプリカ交換の仕方を隣接温度パラメータ間のみの交換にすることができる。

図１０は、隣接温度パラメータ間のみによるレプリカ交換の例を示す図である。
図１０では、レプリカ番号＝１のレプリカに最も小さい温度パラメータ、レプリカ番号＝２のレプリカに次に小さい温度パラメータ、レプリカ番号＝３のレプリカに２番目に大きい温度パラメータ、レプリカ番号＝４のレプリカに最も大きい温度パラメータが設定されているものとする。

たとえば、レプリカ番号＝２のレプリカの作り出す確率密度分布がレプリカ番号＝１，３のレプリカの作り出す確率密度分布とのみ重なり合いがあるのであれば、レプリカ番号＝２のレプリカはレプリカ番号＝４のレプリカとの交換試行は必ず棄却される。このように複数の温度パラメータの値を選ぶことで隣接する温度パラメータが設定されるレプリカ間のみの交換試行を行えばよいことが分かる。厳密には、レプリカ番号＝２，４のレプリカ間の確率密度分布の重なり合いはゼロではない。しかし、計算は常に有限の桁数で行われ、温度パラメータを適切に設定すれば、交換確率を０として扱ってよい状況を作り出すことができる。

一方で、レプリカ交換法においてよく採用される実装は２種類ある。
１つ目は隣接交換と呼ばれる実装である。これは隣接する温度パラメータが設定されるレプリカの確率密度分布だけが重なり合いをもつ条件を想定する。実装も簡単であり、よく使われる条件であるが、温度パラメータの設定については上記の条件を守るように計算者側が事前に予備計算などで温度パラメータの値を決めておくことになる。

２つ目はランダム交換と呼ばれる方法である。これは乱数を用いて任意の２つのレプリカを選択し、全てのレプリカ同士を交換対象とする方法である。この方法は長時間平均を取れば、任意の２つのレプリカの交換試行がなされる。そのため、あるレプリカのエネルギー空間上の確率密度分布が、隣接する温度パラメータ以外が設定されるレプリカの確率密度分布と無視できない大きさの重なり合いをもつ場合でも、詳細つり合いの条件が満たされる特徴がある。どちらの方法を採用しても最小値求解に目的を限れば、既約であることは保証されるため、原理的に解に到達できないということはない。

しかし、効率性を考えた場合、レプリカ間の詳細つり合いの条件を満たさない場合、不変分布条件は満たされない。そのため、全レプリカ系の作り出すエネルギー空間上の確率密度分布を制御するには不利になる可能性がある。ここでの目的はレプリカ交換に必要なレプリカ数を削減することであるから、レプリカ数を最小限度に抑えるため、隣接する温度パラメータが設定されるレプリカ間におけるレプリカ交換で不変分布条件が守られるように温度パラメータの決定が行われる。

このようにして、最小レプリカ数と複数の温度パラメータの値を決定することができる。より具体的な決定方法について、以下に説明する。
まず、レプリカ数が多めに設定され、図６のようにレプリカ交換なしでエネルギー空間上の確率密度分布が求められる。そして、温度パラメータの各値についての確率密度分布の頂点を与えるエネルギーの値とその確率密度分布の広がり（標準偏差）が求められる。これから、温度パラメータの関数として確率密度分布の頂点を与えるエネルギーを求めることができる。一方で確率密度分布の広がりから重なり合いの程度を求めることができる。そして頂点を与えるエネルギーと重なり合いの程度から、交換確率の大小を考慮しつつ、隣接する温度パラメータの値のみが交換されるような複数の温度パラメータを選択することができる。

図５に示したように、べき乗型の遷移確率を用いた場合に得られるエネルギー（期待値）の分布が高温領域では正規分布に近くなることを利用して、たとえば、以下の式（３４）を満たすように温度パラメータであるβ_ｉを選べばよい。

式（３４）において、左辺の１項目は、レプリカ番号＝ｉのレプリカによるエネルギー空間上の確率密度分布の頂点を与えるエネルギー、左辺の２項目はその確率密度分布の標準偏差に所定の係数ｎを乗じた値である。右辺の１項目は、レプリカ番号＝ｉ＋１のレプリカによる確率密度分布の頂点を与えるエネルギー、右辺の２項目はその確率密度分布の標準偏差にｎを乗じた値である。

ｎは確率密度分布同士の重なりの大きさを示す変数である。交換確率を大きくとるのであれば、重なり合いを大きくする必要があるため、レプリカ数は増える。交換確率を小さくしてもよいのであれば、レプリカ数を減らすことができる。

図１１は、エネルギー空間上の確率密度分布の頂点を与えるエネルギーと確率密度分布の標準偏差の例を示す図である。横軸はエネルギー（Ｅ）、縦軸は確率密度ｎ（Ｅ）を表す。なお、図１１には、図６のｍ＝２のべき乗型の遷移確率を用いた場合の確率密度分布の一部が示されている。

図１１には、式（３４）の左辺の１項目である、レプリカ番号＝ｉのレプリカによるエネルギー空間上の確率密度分布の頂点を与えるエネルギーと、その確率密度分布の標準偏差の例が示されている。さらに、図１１には、右辺の２項目である、レプリカ番号＝ｉ＋１のレプリカによる確率密度分布の頂点を与えるエネルギーと、その確率密度分布の標準偏差の例が示されている。

なお、温度パラメータであるＴの計算では、Ｔが低い順に決められていく。
図１２は、低温領域での確率密度分布の振る舞いの例を示す図である。横軸はエネルギー、縦軸は確率密度を表す。

図１２には低温領域でのエネルギー空間上の確率密度分布が表されている。低温領域では温度パラメータの値が低すぎると（図１２の“Ｔ：小”参照）、遷移自体が起こらなくなる。そのため、低エネルギー側の探索効率が落ちてしまう。逆に温度パラメータの値が高くなってしまうと（図１２の“Ｔ：大”参照）、逆に遷移しすぎてしまい低温側の探索効率が落ちてしまう。そのため、一番小さい温度パラメータの値を選択するときには最適な値（たとえば、図１２の“Ｔ：中”）が存在する。この最適な値は、図５に示したように、温度パラメータの関数としてエネルギーを計算したとき、ある温度パラメータの値でエネルギーの期待値が最小値を取ることから決めることができる。

なお、温度パラメータの最小値は予備計算の試行回数を短く取ることから生じるアーチファクトの影響を受ける。しかし、温度パラメータの値を最適化した後の計算においても十分多くの試行回数を取ることは難しいため、見せかけの最小値を与える温度パラメータの値を最小値として取ることにする。

温度パラメータの最小値が決まった後は、式（３４）にしたがって残りの温度パラメータの値を決めていけばよい。このとき、温度パラメータの関数としてのエネルギー関数が補間法などを用いて求められる。標準偏差についても同様に補間法などを用いて求められる。補間法については何を用いてもよいが、低温領域で誤差が大きくなるため、最小二乗法を用いた曲線補間などを用いて平滑化された曲線を求めておけばよい。補間曲線を求めた後は、式（３４）を満たすようなレプリカ番号＝ｉ＋１のレプリカによる確率密度分布の頂点のエネルギーと、その確率密度分布の標準偏差を求めればよい。ここから対応するβ_ｉ＋１、つまりＴ_ｉ＋１を求めることができる。

（ハードウェア構成例）
上記のようなレプリカ交換法や温度パラメータの決定方法については、たとえば以下のようなハードウェアにより実現できる。

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

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

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

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

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

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

次に、情報処理装置２０の機能及び処理手順を説明する。
図１４は、第２の実施の形態の情報処理装置の機能例を示すブロック図である。
情報処理装置２０は、記憶部３０、処理部３１を有する。処理部３１は、制御部３１ａ、設定読込部３１ｂ、スピン初期化部３１ｃ、温度計算部３１ｄ、確率密度計算部３１ｅ、レプリカ交換計算部３１ｆ、結果出力部３１ｇを有する。

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

記憶部３０は、エネルギー情報、スピン情報、レプリカ情報、確率密度情報、問題設定情報、ハミルトニアン情報を記憶する。
エネルギー情報は、計算されたエネルギーの初期値や、これまで計算されたエネルギーの最小値を含む。また、エネルギー情報は、最小値のエネルギーに対応した各状態変数の値を含んでいてもよい。スピン情報は、各状態変数の値を含む。レプリカ情報は、レプリカ交換法を実行するために用いられる情報であり、レプリカ数（Ｎ_{ｒｅｐｌｉｃａ}）、レプリカ交換頻度（Ｎ_ｅｘ）、最低温度を表す温度パラメータの値（Ｔ_ｍｉｎ）、最高温度を表す温度パラメータの値（Ｔ_ｍａｘ）を含む。確率密度情報は、式（２８）の交換確率を計算するための、確率密度の情報（ｎ（β_ｉ，Ｅ_ｊ）など）を含む。確率密度情報は、さらに、たとえば、後述のように確率密度分布をヒストグラムで評価する場合のヒストグラムのビンの数（Ｎ_ｂｉｎ）、確率密度を更新する頻度（Ｎ_ｐｒｏｂ）を含む。問題設定情報は、使用する遷移確率の情報（前述のべき乗型の遷移確率の指数（ｍ）の値）、予備計算のための計算回数（Ｎ_ｐｒｅ）、予備計算後の最適解の求解のための計算回数（Ｎ_ｉｔｅｒ）、スピン初期化法（状態変数の初期値の決め方）の情報を含む。ハミルトニアン情報は、たとえば、式（１）に示したエネルギー関数の重み係数（Ｗ_ｉｊ）、バイアス係数（ｂ_ｉ）、定数（Ｃ）などを含み、前述の問題情報の一例である。

制御部３１ａは、処理部３１の各部を制御する。
設定読込部３１ｂは、記憶部３０から上記の各種情報を、制御部３１ａが理解可能な形式で読み込む。

スピン初期化部３１ｃは、スピン（状態変数）の初期化を行う。
温度計算部３１ｄは、各レプリカに設定する温度パラメータを決定する。
確率密度計算部３１ｅは、式（２８）の交換確率を計算するための、確率密度（ｎ（β_ｉ，Ｅ_ｊ）など）を計算する。

レプリカ交換計算部３１ｆは、レプリカ交換法による求解処理（以下、レプリカ交換処理という）を実行する。
結果出力部３１ｇは、レプリカ交換処理の結果（探索結果）を出力する。結果出力部３１ｇは、たとえば、レプリカ交換処理が所定の終了条件を満たすときに、その時点までに得られた最小のエネルギーとそのエネルギーを与える各状態変数の値を、探索結果として出力する。

図１５は、第２の実施の形態の情報処理装置の一例の処理の流れを示すフローチャートである。
処理が開始すると、まず、設定読込部３１ｂが、記憶部３０から上記の各種情報を、制御部３１ａが理解可能な形式で読み込む（ステップＳ１０）。その後、スピン初期化部３１ｃは、状態変数の初期化を行う（ステップＳ１１）。また、温度計算部３１ｄと確率密度計算部３１ｅによる予備計算が行われる（ステップＳ１２）。ステップＳ１２の処理では、温度計算部３１ｄにより、温度パラメータの計算が行われ、確率密度計算部３１ｅにより、交換確率の計算に用いる確率密度が計算される。計算された確率密度の情報は、記憶部３０に記憶される。

その後、制御部３１ａは、複数のレプリカにそれぞれ異なる温度パラメータの値を設定する（ステップＳ１３）。
そして、レプリカ交換計算部３１ｆによるレプリカ交換処理が行われ（ステップＳ１４）、その結果が出力される（ステップＳ１５）。結果出力部３１ｇは、たとえば、レプリカ交換処理が所定の終了条件（たとえば、計算回数がＮ_ｉｔｅｒに達したこと）を満たすときに、その時点までに得られた最小のエネルギーとそのエネルギーを与える各状態変数の値を、レプリカ交換処理の結果として出力する。

その後、情報処理装置２０の処理が終了する。なお、以下、ステップＳ１４の処理を、ステップＳ１２の予備計算と対比させて本計算と呼ぶ場合もある。
（情報読込処理の例）
図１６は、情報読込処理の一例の処理の流れを示すフローチャートである。

設定読込部３１ｂは、記憶部３０からハミルトニアン情報（式（１）に示したエネルギー関数の重み係数（Ｗ_ｉｊ）、バイアス係数（ｂ_ｉ）、定数（Ｃ））を読み込む（ステップＳ２０）。

また、設定読込部３１ｂは、記憶部３０からＴ_ｍｉｎ、Ｔ_ｍａｘを読み込む（ステップＳ２１）。
さらに、設定読込部３１ｂは、記憶部３０からＮ_{ｒｅｐｌｉｃａ}、Ｎ_ｐｒｅ、Ｎ_ｉｔｅｒ、Ｎ_ｅｘ、Ｎ_ｂｉｎ、Ｎ_ｐｒｏｂを読み込む（ステップＳ２２）。

また、設定読込部３１ｂは、記憶部３０からスピン初期化法を読み込み（ステップＳ２３）、情報読込処理を終える。
なお、図１６に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。

（スピン初期化処理の例）
図１７は、スピン初期化処理の一例の処理の流れを示すフローチャートである。
スピン初期化部３１ｃは、スピン初期化法が指定モードであるか否かを判定する（ステップＳ３０）。スピン初期化部３１ｃは、スピン初期化法が指定モードであると判定した場合、情報処理装置２０の外部から指定された各状態変数の初期値で、全状態変数を初期化し（ステップＳ３１）、スピン初期化処理を終える。

スピン初期化部３１ｃは、スピン初期化法が指定モードではないと判定した場合、スピン初期化法が０モードであるか否かを判定する（ステップＳ３２）。スピン初期化部３１ｃは、スピン初期化法が０モードであると判定した場合、全状態変数を０で初期化し（ステップＳ３３）、スピン初期化処理を終える。スピン初期化部３１ｃは、スピン初期化法が０モードではないと判定した場合、全状態変数を１で初期化し（ステップＳ３４）、スピン初期化処理を終える。

なお、図１７に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
（温度パラメータ計算処理の例）
次に、図１５の予備計算の１つ目の処理である温度パラメータの計算処理の例を説明する。

図１８は、温度パラメータ計算処理の一例の処理の流れを示すフローチャートである。
温度計算部３１ｄは、複数の温度パラメータ（Ｔ）のそれぞれの値について得られるエネルギー空間上の確率密度分布の頂点のエネルギーから、Ｔの関数であるエネルギーの補間曲線（Ｅ_ａｖｅ（Ｔ））を算出する（ステップＳ４０）。複数の温度パラメータ（Ｔ）のそれぞれの値について得られるエネルギー空間上の確率密度分布の頂点のエネルギーは、たとえば図４に示したようなべき乗型の遷移確率を用いた場合のエネルギーの温度依存性のサンプリング結果から得られる。

また、温度計算部３１ｄは、複数の温度パラメータ（Ｔ）のそれぞれの値について得られるエネルギー空間上の確率密度分布の標準偏差から、Ｔの関数である標準偏差の補間曲線（σ（Ｔ））を算出する（ステップＳ４１）。複数の温度パラメータ（Ｔ）のそれぞれの値について得られるエネルギー空間上の確率密度分布の標準偏差は、たとえば図４に示したようなべき乗型の遷移確率を用いた場合のエネルギーの温度依存性のサンプリング結果から得られる。

そして、温度計算部３１ｄは、Ｅ_ａｖｅ（Ｔ）とσ（Ｔ）に基づいて、たとえば、２分法を用いて、式（３４）の関係を満たすように、Ｔ_ｍｉｎからＴ_ｍａｘまでの温度パラメータ（Ｔ）の値を、Ｎ_{ｒｅｐｌｉｃａ}数分決定する（ステップＳ４２）。その後、温度計算部３１ｄは、温度パラメータの計算処理を終える。

なお、ステップＳ４０，Ｓ４１の処理順序は逆であってもよい。
ところで、温度パラメータの計算処理にあまり手間をかけないようにする場合は、前述の隣接交換の条件は厳密に守られなくてもよい。なぜならば、最小値求解問題では平衡状態の存在と、既約性だけが保証されていればよいからである。この条件を厳密に守ることが望ましいのは、統計力学的なシミュレーションで何らかの物理量を計算し、詳細つり合いの条件を厳密に守り、特定のアンサンブルを作り出すことが望ましい場合である。

（確率密度計算処理の一例）
次に、図１５の予備計算の２つ目の処理である確率密度の計算処理の例を説明する。
確率密度（ｎ（β_ｉ，Ｅ_ｊ）など）は、式（２９）のレプリカ交換の交換確率（Ｐ_ｅｘ）を計算するために算出される。確率密度は、温度パラメータの各値を用いたべき乗型の遷移確率を用いたＭＣＭＣ計算の実行時に、温度パラメータの各値について独立なサンプリング計算をすることで容易に計算できる。たとえば、確率密度を求める簡単な方法として、ヒストグラムを用いた近似計算がある。比較的短時間の予備計算を全ての温度パラメータのそれぞれの値について行うと、レプリカ交換なしでのエネルギーの最小値（Ｅ_ｍｉｎ）と最大値（Ｅ_ｍａｘ）を求めることができる。ヒストグラムのビンの数をＮ_ｂｉｎとすると、Ｎ_ｂｉｎ＋１個の点がある。このとき、各ビンの幅を一定とすれば、ｉ番目の点は、以下の式（３５）で表せる。

式（３５）において、ｉ＝０，１，２，…，Ｎ_ｂｉｎである。Ｅ_ｍｉｎ＝Ｅ_{ｂｉｎ，０}、Ｅ_ｍａｘ＝Ｅ_{ｂｉｎ，Ｎ}とすると、区間［Ｅ_{ｂｉｎ，ｉ}，Ｅ_{ｂｉｎ，ｉ＋１}］を定めることができる。

図１９は、確率密度の計算処理の一例の処理の流れを示すフローチャートである。
確率密度計算部３１ｅは、上記の方法で、温度パラメータの各値についてのＥ_ｍｉｎとＥ_ｍａｘを決定し（ステップＳ５０）、式（３５）によりＥ_{ｂｉｎ，ｉ}を決定する（ステップＳ５１）。

そして、確率密度計算部３１ｅは、変数ｋを１にして（ステップＳ５２）、変数ｊを１にして（ステップＳ５３）、変数ｉを０にする（ステップＳ５４）。
その後、確率密度計算部３１ｅは、Ｅ_{ｂｉｎ，ｉ}≦Ｅ_ｊ≦Ｅ_{ｂｉｎ，ｉ＋１}であるか否かを判定する（ステップＳ５５）。Ｅ_ｊは、レプリカ番号＝ｋのレプリカに設定される温度パラメータの値についてのサンプリング計算で得られたＮ_ｄａｔａ個のデータ（エネルギーの値）のうち、ｊ番目のデータである。

確率密度計算部３１ｅは、Ｅ_{ｂｉｎ，ｉ}≦Ｅ_ｊ≦Ｅ_{ｂｉｎ，ｉ＋１}であると判定した場合、レプリカ番号＝ｋのレプリカに設定される温度パラメータの値についての、区間［Ｅ_{ｂｉｎ，ｉ}，Ｅ_{ｂｉｎ，ｉ＋１}］におけるデータ点数を示すｎ_ｉ ^ｋを＋１する（ステップＳ５６）。

確率密度計算部３１ｅは、ステップＳ５６の後、またはＥ_{ｂｉｎ，ｉ}≦Ｅ_ｊ≦Ｅ_{ｂｉｎ，ｉ＋１}ではないと判定した場合、ｉ＝Ｎ_ｂｉｎであるか否かを判定する（ステップＳ５７）。確率密度計算部３１ｅは、ｉ＝Ｎ_ｂｉｎではないと判定した場合、ｉを＋１して（ステップＳ５８）、ステップＳ５５からの処理を繰り返す。

確率密度計算部３１ｅは、ｉ＝Ｎ_ｂｉｎであると判定した場合、ｊ＝Ｎ_ｄａｔａであるか否かを判定する（ステップＳ５９）。確率密度計算部３１ｅは、ｊ＝Ｎ_ｄａｔａではないと判定した場合、ｊを＋１して（ステップＳ６０）、ステップＳ５４からの処理を繰り返す。

確率密度計算部３１ｅは、ｊ＝Ｎ_ｄａｔａであると判定した場合、ｋ＝Ｎ_{ｒｅｐｌｉｃａ}であるか否かを判定する（ステップＳ６１）。確率密度計算部３１ｅは、ｋ＝Ｎ_{ｒｅｐｌｉｃａ}ではないと判定した場合、ｋを＋１して（ステップＳ６２）、ステップＳ５３からの処理を繰り返す。

確率密度計算部３１ｅは、ｋ＝Ｎ_{ｒｅｐｌｉｃａ}であると判定した場合、ｋ＝１（ステップＳ６３）、ｉ＝０とする（ステップＳ６４）。
その後、確率密度計算部３１ｅは、ステップＳ６２までの処理で得たｎ_ｉ ^ｋを、ｎ_ｉ ^ｋ／Ｎ_ｄａｔａで更新する（ステップＳ６５）。

そして、確率密度計算部３１ｅは、ｉ＝Ｎ_ｂｉｎであるか否かを判定する（ステップＳ６６）。確率密度計算部３１ｅは、ｉ＝Ｎ_ｂｉｎではないと判定した場合、ｉを＋１して（ステップＳ６７）、ステップＳ６５からの処理を繰り返す。

確率密度計算部３１ｅは、ｉ＝Ｎ_ｂｉｎであると判定した場合、ｋ＝Ｎ_{ｒｅｐｌｉｃａ}であるか否かを判定する（ステップＳ６８）。確率密度計算部３１ｅは、ｋ＝Ｎ_{ｒｅｐｌｉｃａ}ではないと判定した場合、ｋを＋１して（ステップＳ６９）、ステップＳ６４からの処理を繰り返す。

確率密度計算部３１ｅは、ｋ＝Ｎ_{ｒｅｐｌｉｃａ}であると判定した場合、確率密度の計算処理を終える。
なお、図１９に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。

上記の処理によって得られたｎ_ｉ ^ｋから、レプリカ交換の交換確率（Ｐ_ｅｘ）を計算するために用いられる確率密度が得られる。
なお、上記の温度パラメータの各値についてのＥ_ｍｉｎとＥ_ｍａｘは本計算時に更新され、上記の確率密度は更新される。

（レプリカ交換処理（本計算）の例）
本計算では、ステップＳ１２の処理で決定された複数の温度パラメータの値の何れかが設定された複数のレプリカのそれぞれにおいて、べき乗型の遷移確率を用いたＭＣＭＣ計算が行われる。

ＭＣＭＣ計算の途中で前述のＥ_ｍｉｎまたはＥ_ｍａｘが更新された場合、確率密度計算部３１ｅは、更新されたＥ_ｍｉｎまたはＥ_ｍａｘを用いて前述のヒストグラムを更新することで、交換確率の計算に用いられる確率密度を更新する。

図２０は、確率密度の更新処理の一例の流れを示すフローチャートである。
確率密度計算部３１ｅは、ＭＯＤ（Ｎ_ｓｔｅｐ，Ｎ_ｐｒｏｂ）＝０であるか否かを判定する（ステップＳ７０）。Ｎ_ｓｔｅｐは本計算の現在の計算回数（ステップ数）であり、Ｎ_ｐｒｏｂは、確率密度の更新頻度を示すステップ数である。ステップＳ７０の処理では、Ｎ_ｓｔｅｐがＮ_ｐｒｏｂの倍数であるか否かが判定される。

ヒストグラムの更新頻度は少なくてよい。Ｎ_ｐｒｏｂは、少なくともサンプリング頻度を示す計算回数よりも十分大きく設定されることが望ましい。たとえば、サンプリングが、本計算の計算回数の１０００回に１回行われるなら、ヒストグラムの更新頻度は１０００回のサンプリングごとに１回行われるようにＮ_ｐｒｏｂが設定される。

確率密度計算部３１ｅは、ＭＯＤ（Ｎ_ｓｔｅｐ，Ｎ_ｐｒｏｂ）＝０ではないと判定した場合、後述の処理によりヒストグラムの最小値（Ｅ_ｍｉｎ）または最大値（Ｅ_ｍａｘ）の更新を行い（ステップＳ７１）、処理を終える。

確率密度計算部３１ｅは、ＭＯＤ（Ｎ_ｓｔｅｐ，Ｎ_ｐｒｏｂ）＝０であると判定した場合、ヒストグラムの更新を行い（ステップＳ７２）、処理を終える。ステップＳ７２の処理では、ヒストグラムを更新するタイミングにおける最新のＥ_ｍｉｎまたはＥ_ｍａｘが用いられる。なお、確率密度計算部３１ｅは、Ｅ_ｍｉｎが更新された場合、ヒストグラムにおいて区間［Ｅ_{ｂｉｎ，０}，Ｅ_{ｂｉｎ，１}］だけ更新し、Ｅ_ｍａｘが更新された場合、ヒストグラムにおいて区間［Ｅ_{ｂｉｎ，Ｎ－１}，Ｅ_{ｂｉｎ，Ｎ}］だけ更新する。これにより、ヒストグラム全体を更新する場合よりも計算量を抑えられる。

図２１は、Ｅ_ｍｉｎとＥ_ｍａｘの更新処理の一例の流れを示すフローチャートである。
確率密度計算部３１ｅは、各レプリカにおけるＭＣＭＣ計算により状態遷移が生じるたびに以下の処理を行う。

確率密度計算部３１ｅは、ＭＣＭＣ計算の繰り返し計算の途中で得られる現在の状態変数の値に対応したエネルギー（Ｅ_ｎｏｗ）が、Ｅ_ｎｏｗ＜Ｅ_ｍｉｎであるか否かを判定する（ステップＳ８０）。

確率密度計算部３１ｅは、Ｅ_ｎｏｗ＜Ｅ_ｍｉｎであると判定した場合、Ｅ_ｍｉｎ＝Ｅ_ｎｏｗに更新する（ステップＳ８１）。確率密度計算部３１ｅは、ステップＳ８１の処理後、またはＥ_ｎｏｗ＜Ｅ_ｍｉｎではないと判定した場合、Ｅ_ｎｏｗ＞Ｅ_ｍａｘであるか否かを判定する（ステップＳ８２）。

確率密度計算部３１ｅは、Ｅ_ｎｏｗ＞Ｅ_ｍａｘであると判定した場合、Ｅ_ｍａｘ＝Ｅ_ｎｏｗに更新する（ステップＳ８３）。確率密度計算部３１ｅは、ステップＳ８３の処理後、またはＥ_ｎｏｗ＞Ｅ_ｍａｘではないと判定した場合、Ｅ_ｍｉｎとＥ_ｍａｘの更新処理を終える。

なお、図２１に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
図２０、図２１のような確率密度の更新処理を行うことで、本計算の現在の計算回数が増えていくと、確率密度のサンプリング精度も上がっていく。

本計算では、レプリカ交換計算部３１ｆは、レプリカ交換頻度を示すＮ_ｅｘごとに、式（２９）の交換確率に基づいて、レプリカ間における温度パラメータの値の交換を行う。なお、温度パラメータの値の交換の代わりに、状態（Ｎ個の状態変数の値）を交換してもよい。

図２２は、レプリカ交換処理の一例の流れを示すフローチャートである。
レプリカ交換計算部３１ｆは、ＭＯＤ（Ｎ_ｓｔｅｐ，Ｎ_ｅｘ）＝０であるか否かを判定する（ステップＳ９０）。ステップＳ９０の処理では、現在の計算回数であるＮ_ｓｔｅｐがレプリカ交換頻度であるＮ_ｅｘの倍数であるか否かが判定される。

レプリカ交換計算部３１ｆは、ＭＯＤ（Ｎ_ｓｔｅｐ，Ｎ_ｅｘ）＝０であると判定した場合、ＭＯＤ（Ｎ_{ｔｏｔ＿ｅｘ}，２）＝０であるか否かを判定する（ステップＳ９１）。レプリカ交換計算部３１ｆは、ＭＯＤ（Ｎ_ｓｔｅｐ，Ｎ_ｅｘ）＝０ではないと判定した場合、レプリカ交換処理を終える。ステップＳ９１の処理では、現在のレプリカ交換回数を示すＮ_{ｔｏｔ＿ｅｘ}が偶数であるか否かが判定される。

レプリカ交換計算部３１ｆは、ＭＯＤ（Ｎ_{ｔｏｔ＿ｅｘ}，２）＝０であると判定した場合、Ｔ_ｏｄｄとＴ_{ｏｄｄ＋１}が設定されるレプリカのペアを交換候補として選択する（ステップＳ９２）。Ｔ_ｏｄｄは、ステップＳ１２の処理で計算された温度パラメータ（Ｔ）の値を、小さい順に並べたときに、奇数番目の温度パラメータの値を示す。たとえば、Ｔ_１，Ｔ_２，Ｔ_３，Ｔ_４，Ｔ_５，…の順に温度パラメータの値が配列されているとする。この場合、Ｔ_１が設定されているレプリカとＴ_２が設定されているレプリカのペア、Ｔ_３が設定されているレプリカとＴ_４が設定されているレプリカのペアが交換候補に含まれる。

レプリカ交換計算部３１ｆは、ＭＯＤ（Ｎ_{ｔｏｔ＿ｅｘ}，２）＝０ではないと判定した場合、Ｔ_ｅｖｅｎとＴ_{ｅｖｅｎ＋１}が設定されるレプリカのペアを交換候補として選択する（ステップＳ９３）。Ｔ_ｅｖｅｎは、温度パラメータ（Ｔ）の値を、小さい順に並べたときに、偶数番目の温度パラメータの値を示す。たとえば、Ｔ_１，Ｔ_２，Ｔ_３，Ｔ_４，Ｔ_５，…の順に温度パラメータの値が配列されているとする。この場合、Ｔ_２が設定されているレプリカとＴ_３が設定されているレプリカのペア、Ｔ_４が設定されているレプリカとＴ_５が設定されているレプリカのペアが交換候補に含まれる。

次に、レプリカ交換計算部３１ｆは、交換候補のペアを１つ選択し（ステップＳ９４）、区間［０，１］の値をもつ乱数Ｒを発生させる（ステップＳ９５）。そして、レプリカ交換計算部３１ｆは、式（２９）の交換確率であるＰ_ｅｘが、Ｐ_ｅｘ≧Ｒであるか否かを判定する（ステップＳ９６）。

レプリカ交換計算部３１ｆは、Ｐ_ｅｘ≧Ｒであると判定した場合、選択したペアのレプリカ間で設定されている温度パラメータの値を交換することでレプリカ交換を実行する（ステップＳ９７）。

レプリカ交換計算部３１ｆは、ステップＳ９７の処理後、または、Ｐ_ｅｘ≧Ｒではないと判定した場合、ステップＳ９２またはステップＳ９３の処理で選択された全交換候補が、ステップＳ９４の処理で全て選択したか否かを判定する（ステップＳ９８）。

レプリカ交換計算部３１ｆは、全交換候補を選択していないと判定した場合、ステップＳ９４からの処理を繰り返し、全交換候補を選択したと判定した場合、１回のレプリカ交換処理を終える。

なお、図２２に示した処理の流れは一例であって、適宜処理の順序が入れ替えられていてもよい。
また、図２２に示した処理は、隣接交換によるレプリカ交換処理であるが、ランダム交換が行われる場合、乱数で交換対象のレプリカのペアが選択されるように変更すればよい。ペアが決定された後の計算手続きは隣接交換の場合と同じである。

ここで、理論面の補足をしておく。式（２６）の任意の２つのレプリカ間の詳細つり合い条件は不変分布条件を課している。
べき乗型の遷移確率が用いられる場合、最小値の求解が目的であるため、不変分布条件は理論的に必須ではない。しかし、べき乗型の遷移確率が何らかの分布を作り出し、その分布が不変分布になるような拘束条件を付けることはサンプリング空間を一定に保つために望ましい。なぜなら、予備計算では、温度パラメータの各値における確率密度分布が見積もられ、その見積もりにしたがって、なるべく広くサンプリングができるように複数の温度パラメータの値とレプリカ数が定義される。これは見積もりに用いられる温度パラメータの各値が設定されるレプリカにおいて得られる確率密度分布がレプリカ交換によって変わらないことを暗に仮定している。そのため、最小値求解問題において、べき乗型の遷移確率を用いた場合においても、不変分布条件を課すことに合理性はある。

（効果）
図２３は、温度パラメータの値の交換の様子を示す図である。横軸は計算回数を表し、縦軸は温度パラメータ（Ｔ）を表す。

図２３では、ボルツマン分布に基づく遷移確率を用いてレプリカ交換処理が行われた場合の温度パラメータの値の交換の様子と、べき乗型の遷移確率（ｍ＝３）を用いてレプリカ交換処理が行われた場合の温度パラメータの値の交換の様子が示されている。

なお、ボルツマン分布に基づく遷移確率を用いたレプリカ交換処理における複数の温度パラメータの値は、最適化されたとみなせるものが用いられている。これに対して、べき乗型の遷移確率を用いたレプリカ交換処理における複数の温度パラメータの値については最適化がされていない。その理由は、効果を検証するために、べき乗型の遷移確率を用いたレプリカ交換処理を、ボルツマン分布に基づく遷移確率を用いたレプリカ交換処理よりも不利な条件下に置いたからである。

べき乗型の遷移確率を用いたレプリカ交換処理における複数の温度パラメータの値について、最低温度がＴ＝５になっており、ボルツマン分布に基づく遷移確率を用いた場合よりも高い温度になっているのはｍの値を比較的大きいｍ＝３にとっているからである。Ｔ＝５は、ｍ＝３の系では十分低温とみなせる値である（ｍ＝１．００１の場合はＴ＝０．１程度が十分低温とみなせる値になる）。

このようにｍの値に応じて、十分低温とみなせる温度パラメータの値が、最低温度を示すＴ_ｍｉｎ、十分高温とみなせる温度パラメータの値が、最高温度を示すＴ_ｍａｘとして決定される。

図２３では、レプリカ交換された隣接温度間に線が表示されている。べき乗型の遷移確率を用いた場合、線が密になっており、全ての温度帯でレプリカ交換が実行されている様子がわかる。一方、ボルツマン分布に基づく遷移確率を用いた場合、高温領域になるほど空白が多くなっており、レプリカ交換が実行されにくくなっている様子がわかる。

図２４は、レプリカに設定される温度パラメータの値の変化を示す図である。横軸は計算回数（ステップ数）を表し、縦軸は温度パラメータ（Ｔ）を表す。
図２４では、ボルツマン分布に基づく遷移確率を用いてレプリカ交換処理が行われた場合のレプリカ番号＝０，２５のレプリカに設定される温度パラメータの値の変化が示されている。また、べき乗型の遷移確率（ｍ＝３）を用いてレプリカ交換処理が行われた場合のレプリカ番号＝４，２５のレプリカに設定される温度パラメータの値の変化が示されている。複数の温度パラメータの値の設定に関しては、図２３と同様である。

図２４のように、ボルツマン分布に基づく遷移確率を用いた場合、レプリカ番号＝０のレプリカには、比較的小さい温度パラメータの値が設定され、レプリカ番号＝２５のレプリカには、比較的大きい温度パラメータの値が設定されている。すなわち、低温領域と高温領域が分離していることがわかる。これは相転移点近傍の温度パラメータの値の調節が難しいことによる。交換確率はエネルギー差に関して指数関数的に小さくなる。相転移点近傍ではエネルギー差が急激に変化し、交換しづらい傾向がある。そのため、相転移点近傍を効率的に超えることが難しくなっているのである。

一方、べき乗型の遷移確率を用いた場合、複数の温度パラメータの値を最適化していないにも関わらず、レプリカは低温領域と高温領域を行き来できることが分かる。このことは、図５の温度パラメータ（Ｔ）の関数としてのエネルギーの図から理解できる。べき乗型の遷移確率を用いた場合、ボルツマン分布に基づく遷移確率を用いた場合と比較すると、Ｔの関数としてエネルギーは急激に増加していない。つまり、レプリカ交換しやすくなっているのである。これはべき乗型の遷移確率を用いたレプリカ交換が温度パラメータの設定に対してロバストであることを意味する。

このようにして定性的にはべき乗型の遷移確率を用いたレプリカ交換処理が有用であることが分かるが、定量的な性能改善を論じるために、定量評価指標としてトンネル時間を採用する。トンネル時間とは１つのレプリカがＴ_ｍｉｎからＴ_ｍａｘを経由して、再度Ｔ_ｍｉｎまで戻る時間である。逆にＴ_ｍａｘからＴ_ｍｉｎを経由して再度Ｔ_ｍａｘに戻る時間をトンネル時間としてもよい。

この指標が定量指標になるのは、レプリカ交換法のアルゴリズムが導入されたモチベーションが、高温領域を経由させてエネルギーランドスケープ上の大域的構造を効率的に変え、サンプリング効率を上げることにあることによる。そのため、効率的にＴ_ｍａｘとＴ_ｍｉｎを行き来できればそれだけ効率的と考えられるためである。

図２５は、ボルツマン分布に基づく遷移確率を用いた場合とべき乗型の遷移確率を用いた場合のレプリカ交換処理時のトンネル時間の比較結果の例を示す図である。横軸はトンネル時間を計算回数（ステップ数）で表したものであり、縦軸は頻度を表す。なお、べき乗型の遷移確率は、ｍ＝３のものを使用した。

ボルツマン分布に基づく遷移確率を用いたレプリカ交換処理の場合、トンネル時間の平均値はおおよそ、１５００００ステップであった。これに対して、べき乗型の遷移確率を用いたレプリカ交換処理の場合、トンネル時間の平均値はおおよそ９２０００ステップであった。比を取ると１．６３倍程度であるため、トンネル時間ベースの比較であれば約６３％の性能向上になっていた。

ただし、ボルツマン分布に基づく遷移確率を用いた場合も繰り返し最適化を行い、労力を考慮しなければトンネル時間を短くすることは可能である。べき乗型の遷移確率を採用する最大のメリットはレプリカ交換処理時に温度パラメータの値に対してロバストになるため、温度パラメータの各値を決定する労力を削減できる点にある。つまり、手抜きをしてもある程度の性能が得られることに意義がある。最適な温度パラメータの値を問題ごとに毎回決定することは大変であるからである。結果を得るための計算機における最終的な計算時間が１０分の１になったとしても、最適な温度パラメータの値を得るための計算や準備に１０倍の時間がかかってしまえば全体最適化にはならない。上記の手法はべき乗型の遷移確率のＴ_ｍｉｎとＴ_ｍａｘを決定してしまえば、温度パラメータの値に対してロバストな方法であるため、いい加減な設定をしても比較的よい性能が得られる。なぜならば、レプリカ交換の効率を落とす要因を意図的に排除したからである。

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

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

１０ 情報処理装置
１１ 記憶部
１２ 処理部

Claims (9)

1. 問題を変換した評価関数の情報を取得し、
   レプリカ交換法による最適解の求解処理に用いる互いに異なる複数の温度パラメータの値を決定し、
   前記複数の温度パラメータの値をそれぞれ複数のレプリカの何れかに１つずつ設定し、
   前記評価関数に含まれる複数の状態変数の何れかの値が変化することによる前記評価関数の値の変化分と前記複数の温度パラメータの何れかの値とに基づいて得られる第１の遷移確率であって、温度パラメータの値の変化に対する前記評価関数の値の変化がボルツマン分布に基づく第２の遷移確率を用いた場合よりも緩やかになる前記第１の遷移確率にしたがって、前記複数の状態変数の何れかの値の更新を繰り返す更新処理を、前記複数のレプリカのそれぞれについて互いに独立に行うとともに、前記第１の遷移確率によって得られる確率分布の不変分布条件を満たす交換確率にしたがって、前記複数のレプリカの間で、前記複数のレプリカのそれぞれに設定された前記複数の温度パラメータの何れかの値、または前記複数のレプリカのそれぞれにおける前記複数の状態変数の値を交換する交換処理を繰り返すことで、前記求解処理を実行する、
   処理をコンピュータに実行させる最適化プログラム。
2. 前記第１の遷移確率は、前記変化分と前記複数の温度パラメータの何れかの値との積に１を加えた値のｍ（ｍ＞１）乗の逆数で表される、請求項１に記載の最適化プログラム。
3. 前記ｍは４以下である、請求項２に記載の最適化プログラム。
4. 前記複数の温度パラメータの値のうち、第１の値に基づいて得られる前記第１の遷移確率を用いた前記更新処理により得られる、前記評価関数の値の第１の確率密度分布の頂点を与える前記評価関数の値と、前記第１の確率密度分布の標準偏差に所定の係数を乗じた値との和が、前記複数の温度パラメータの値のうち、第２の値に基づいて得られる前記第１の遷移確率を用いた前記更新処理により得られる、前記評価関数の値の第２の確率密度分布の頂点を与える前記評価関数の値から、前記第２の確率密度分布の標準偏差に前記係数を乗じた値を引いた値に等しくなるように、前記第１の値と前記第２の値を決定する、処理を前記コンピュータに実行させる請求項１乃至３の何れか一項に記載の最適化プログラム。
5. 前記複数のレプリカのうち、前記複数の温度パラメータの値の１つである第３の値が設定される第１のレプリカと、前記複数の温度パラメータの値の１つである第４の値が設定される第２のレプリカとの間で前記交換処理が行われるときの前記交換確率は、前記交換処理の前後の、前記第１のレプリカにおける前記評価関数の値の第１の確率密度と前記第２のレプリカにおける前記評価関数の値の第２の確率密度の積の比で表される、請求項１乃至４の何れか一項に記載の最適化プログラム。
6. 前記第１の確率密度または前記第２の確率密度は、ヒストグラムを用いた近似計算により計算される、請求項５に記載の最適化プログラム。
7. 前記第１の確率密度または前記第２の確率密度は、前記第１のレプリカまたは前記第２のレプリカにおける前記評価関数の値の最小値または最大値が前記求解処理において更新された場合、前記最小値または前記最大値の更新後の値に基づいて更新される、請求項５または６の何れか一項に記載の最適化プログラム。
8. コンピュータが、
   問題を変換した評価関数の情報を取得し、
   レプリカ交換法による最適解の求解処理に用いる互いに異なる複数の温度パラメータの値を決定し、
   前記複数の温度パラメータの値をそれぞれ複数のレプリカの何れかに１つずつ設定し、
   前記評価関数に含まれる複数の状態変数の何れかの値が変化することによる前記評価関数の値の変化分と前記複数の温度パラメータの何れかの値とに基づいて得られる第１の遷移確率であって、温度パラメータの値の変化に対する前記評価関数の値の変化がボルツマン分布に基づく第２の遷移確率を用いた場合よりも緩やかになる前記第１の遷移確率にしたがって、前記複数の状態変数の何れかの値の更新を繰り返す更新処理を、前記複数のレプリカのそれぞれについて互いに独立に行うとともに、前記第１の遷移確率によって得られる確率分布の不変分布条件を満たす交換確率にしたがって、前記複数のレプリカの間で、前記複数のレプリカのそれぞれに設定された前記複数の温度パラメータの何れかの値、または前記複数のレプリカのそれぞれにおける前記複数の状態変数の値を交換する交換処理を繰り返すことで、前記求解処理を実行する、
   最適化方法。
9. 問題を変換した評価関数の情報を記憶する記憶部と、
   前記情報を取得し、レプリカ交換法による最適解の求解処理に用いる互いに異なる複数の温度パラメータの値を決定し、前記複数の温度パラメータの値をそれぞれ複数のレプリカの何れかに１つずつ設定し、前記評価関数に含まれる複数の状態変数の何れかの値が変化することによる前記評価関数の値の変化分と前記複数の温度パラメータの何れかの値とに基づいて得られる第１の遷移確率であって、温度パラメータの値の変化に対する前記評価関数の値の変化がボルツマン分布に基づく第２の遷移確率を用いた場合よりも緩やかになる前記第１の遷移確率にしたがって、前記複数の状態変数の何れかの値の更新を繰り返す更新処理を、前記複数のレプリカのそれぞれについて互いに独立に行うとともに、前記第１の遷移確率によって得られる確率分布の不変分布条件を満たす交換確率にしたがって、前記複数のレプリカの間で、前記複数のレプリカのそれぞれに設定された前記複数の温度パラメータの何れかの値、または前記複数のレプリカのそれぞれにおける前記複数の状態変数の値を交換する交換処理を繰り返すことで、前記求解処理を実行する処理部と、
   を有する情報処理装置。

Images