JP2021043503A - 最適化装置、最適化方法及び最適化プログラム - Google Patents

Abstract

【課題】レプリカ交換確率の均一性を高める。【解決手段】最適化装置１０は、温度変化に応じたエネルギーの変化度合いを示す比熱と温度との関係を示す比熱情報を取得する。そして最適化装置１０は、比熱情報に基づいて、レプリカ交換法でエネルギー関数の最適値を求解する際に複数のレプリカそれぞれに設定する複数の温度値を決定する。【選択図】図１

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

記憶部１１は、エネルギー関数の最小値を求める最小値求解問題を示す求解問題情報１１ａと熱平衡状態情報１１ｂとが含まれる。求解問題情報１１ａには、たとえば求解対象のエネルギー関数やエネルギー関数に設定する状態変数の初期値などが含まれる。熱平衡状態情報１１ｂには、たとえば複数の標本温度値（Ｔ₁₁，Ｔ₁₂，・・・）に対応付けられた、標本温度値（Ｔ₁₁，Ｔ₁₂，・・・）ごとの熱平衡状態でのエネルギー値（Ｅ₁，Ｅ₂，・・・）が含まれる。エネルギー値（Ｅ₁，Ｅ₂，・・・）は、たとえばエネルギー関数を複数の標本温度値（Ｔ₁₁，Ｔ₁₂，・・・）それぞれで解くことで得られる。

熱平衡状態情報１１ｂは、１つの標本温度値に対応付けられた、その標本温度値において熱平衡状態で繰り返しエネルギー関数を計算することで算出された複数のエネルギー値を含んでいてもよい。また熱平衡状態情報１１ｂは、１つの標本温度値に対応付けられた、その標本温度値において熱平衡状態で繰り返しエネルギー関数を計算することで算出された複数のエネルギー値の平均値を含んでいてもよい。

処理部１２は、まず予備計算として、複数の標本温度値（Ｔ₁₁，Ｔ₁₂，・・・）それぞれにおいて、エネルギー関数の計算を繰り返し、熱平衡状態になった後のエネルギー値（Ｅ₁，Ｅ₂，・・・）を熱平衡状態情報１１ｂとして記憶部１１に格納する（ステップＳ１）。

次に処理部１２は、温度変化に応じたエネルギーの変化度合いを示す比熱（Ｃ＝ｄＥ／ｄＴ［Ｋ］）と温度（Ｔ）との関係を示す比熱情報を取得する（ステップＳ２）。たとえば処理部１２は、熱平衡状態情報１１ｂに基づいて、比熱（Ｃ＝ｄＥ／ｄＴ［Ｋ］）と温度（Ｔ）との関係を示す比熱情報を算出する。たとえば処理部１２は、複数の標本温度値（Ｔ₁₁，Ｔ₁₂，・・・）を値の大きさで並べたときに隣接する２つの標本温度値（Ｔ₁₁，Ｔ₁₂，・・・）それぞれに対応する比熱の間を補完する補間式を、比熱情報として生成する。

次に処理部１２は、算出した比熱情報に基づいて、レプリカ交換法でエネルギー関数の最小値を求解する際に複数のレプリカそれぞれに設定する複数の温度値（Ｔ₂₁，Ｔ₂₂，Ｔ₂₃，・・・）を決定する（ステップＳ３）。たとえば処理部１２は、第１の温度帯域における温度値同士の間隔が、第１の温度帯域よりも比熱が小さい第２の温度帯域における温度値同士の間隔よりも狭くなるように、複数のレプリカそれぞれに設定する複数の温度値を決定する。

より具体的には、処理部１２は、比熱（Ｃ₁，Ｃ₂，Ｃ₃，・・・）が大きいほど値が小さくなる温度刻み幅算出式に基づいて、温度刻み幅（ΔＴ₁，ΔＴ₂，ΔＴ₃，・・・）を算出する。たとえば処理部１２は、比熱の逆数に定数ｋ倍（ｋは正の実数）した値の平方根を温度値（Ｔ₂₁，Ｔ₂₂，Ｔ₂₃，・・・）に乗算する式を、その温度値（Ｔ₂₁，Ｔ₂₂，Ｔ₂₃，・・・）での温度刻み幅算出式とする。なお定数ｋの算出方法は、第２の実施形態で詳細に説明する（式（１７）等参照）。

たとえば処理部１２は、温度が低い方からｍ番目（ｍは１以上の整数）の温度値における比熱に応じた温度刻み幅を、温度刻み幅算出式に基づいて算出する。そして処理部１２は、ｉ番目の温度値よりも算出した温度刻み幅だけ高い温度を、温度が低い方からｍ＋１番目の温度値に決定する。

たとえば最も低い温度の温度値Ｔ₂₁は、初期値として予め指定されている。そのとき処理部１２は、温度値Ｔ₂₁における温度刻み幅ΔＴ₁＝Ｔ₂₁（ｋ／Ｃ₁）^1/2を計算する。次に処理部１２は、温度値Ｔ₂₁に温度刻み幅ΔＴ₁＝Ｔ₂₁（ｋ／Ｃ₁）^1/2を加算した値（Ｔ₂₁＋ΔＴ₁）を、２番目の温度値Ｔ₂₂に決定する。以後同様に、処理部１２は、温度が低い方から順に、温度値を決定していく。

そして処理部１２は、複数のレプリカそれぞれに複数の温度値（Ｔ₂₁，Ｔ₂₂，Ｔ₂₃，・・・）を設定し、レプリカ交換法によりエネルギー関数の最小値を求解する（ステップＳ４）。

これにより、レプリカ交換確率の均一性を高めることができる。すなわち、レプリカ間の温度間隔が一定の場合、比熱が高い温度帯域では、比熱が低い温度帯域よりもレプリカ交換確率が低くなる。他方、レプリカ間の温度間隔は狭いほど、レプリカ交換確率が高くなる。そこで最適化装置１０では、比熱が高い温度帯域ほどレプリカ間の温度間隔を狭くすることで、レプリカ交換確率の均一性を高めている。

レプリカ交換確率の均一性が向上することで、レプリカは温度空間内を適切にランダムウォークすることができるようになる。その結果、サンプリング空間が広がり、エネルギー関数の最小値の求解において、求解の精度が向上する。

なお、処理部１２は、決定した温度値により、標本温度値の適性化を行うこともできる。たとえば処理部１２は、複数の温度値を決定後、複数の標本温度値を、決定した温度値と同じ値に変更する。そして処理部１２は、変更後の複数の標本温度値それぞれにおいて、エネルギー関数の計算を繰り返し、熱平衡状態になった後のエネルギー値で熱平衡状態情報１１ｂを更新する。これにより、温度値における比熱の誤差を削減することができ、精度の高い温度値を得ることができる。温度値の精度が向上することで、エネルギー関数の最小値の求解も高精度に行うことが可能となる。

〔第２の実施の形態〕
次に第２の実施の形態について説明する。第２の実施の形態は、問題を変換したイジング型のエネルギー関数の最小値（または最小値が得られる状態変数の値の組合せ）を、最適化装置を用いてレプリカ交換法により求解するものである。

図２は、最適化装置のハードウェアの一例を示す図である。最適化装置１００は、たとえばコンピュータであり、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。

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

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

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

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

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

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

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

ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。

最適化装置１００は、以上のようなハードウェア構成によって、第２の実施の形態の処理機能を実現することができる。なお、第１の実施の形態に示した最適化装置１０も、図２に示した最適化装置１００と同様のハードウェアにより実現することができる。

最適化装置１００は、たとえばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態の処理機能を実現する。最適化装置１００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。たとえば、最適化装置１００に実行させるプログラムをストレージ装置１０３に格納しておくことができる。プロセッサ１０１は、ストレージ装置１０３内のプログラムの少なくとも一部をメモリ１０２にロードし、プログラムを実行する。また最適化装置１００に実行させるプログラムを、光ディスク２４、メモリ装置２５、メモリカード２７などの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、たとえばプロセッサ１０１からの制御により、ストレージ装置１０３にインストールされた後、実行可能となる。またプロセッサ１０１が、可搬型記録媒体から直接プログラムを読み出して実行することもできる。

このようなハードウェアの最適化装置１００を用いて、イジング型のエネルギー関数の最小値、および最小値を得る状態変数の組合せを求めることができる。求解対象のエネルギー関数は、ハミルトニアンによって表すことができる。

イジング型のエネルギー関数（Ｈ（｛ｘ｝））は、たとえば、以下の式（１）で定義される。

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

最適化装置１００は、Ｎ個の状態変数のエネルギー関数（Ｈ（｛ｘ｝））の最小値を求める。たとえば最適化装置１００は、マルコフ連鎖モンテカルロ計算を実行し、統計分布としてボルツマン分布の下で実行する。ボルツマン分布を計算に使用するのは、イジングモデルが元々磁性体の物性をシミュレーションするために考案されたモデルであるためである。最適化装置１００は、シミュレーションにおいて、以下のようにして実際の計算を行う。

まず最適化装置１００は、疑似乱数発生器を用いて、反転させる状態変数ｘ_iをランダムに選ぶ。最適化装置１００は、選択後、エネルギー差ΔＥを計算する。選ばれた状態変数をｘ_kとすると、その状態変数ｘ_ｋの値が変化することによるエネルギー関数の値の変化（エネルギー差（ΔＥ））は、以下の式で表される。

式（２）の「１−２ｘ_ｋ」はｘ_ｋの変化分（Δｘ_ｋ）を表す。ｘ_ｋが１から０に変化する場合、１−２ｘ_ｋ＝−１であり、ｘ_ｋが０から１に変化する場合、１−２ｘ_ｋ＝１である。また、ｈ_ｋはローカルフィールドと呼ばれ、以下の式（３）で表せる。

式（２）、式（３）により、エネルギー差ΔＥが求まるので、最適化装置１００は、選択した状態変数ｘ_kを変化させた新たな状態を受け入れるか否かを決定する。たとえば最適化装置１００は、メトロポリス法に従い、式（４）の形の遷移確率を用いて新しい状態を受け入れるかどうか決定する。

式（４）において、π_ｉ（｛ｘ｝）は、｛ｘ｝（各状態変数の値を示す）で規定される状態が実現する確率を示す確率分布である。Ｐ_ｉ→ｊは、状態ｉから状態ｊへ遷移する遷移確率を表す。βは逆温度（温度パラメータの値の逆数）である。

最適化装置１００は、このような確率的な状態遷移を繰り返す。これにより、与えられた温度における熱平衡状態が達成される。たとえば処理部１２が、温度を十分下げた状態でシミュレーションを実行することで、エネルギー最小値を求めることができる。

ＳＡ法では、シミュレーションの進行とともに温度Ｔを下げていき、温度Ｔを０に近づけていくこととなる。温度Ｔを０に近づけることでエネルギーが上がる方向への遷移が抑えられるため、式（１）のエネルギーの停留点が求まる。これは停留点であり、最小値とは限らない。そのため、ＳＡ法では初期条件を変えて多数回の試行を行うことで、ヒューリスティックに最適解を求める戦略が取られる。

なお、式（１）を、ボルツマン分布を用いてメトロポリス法で数値的に解く場合、最適解を求めることについて以下の阻害要因がある。
［阻害要因１］エネルギー関数の局所最適解が無数にあり、容易にトラップされ、抜け出せない。そのため、求めたい解とはかけ離れた解が求まってしまう。
［阻害要因２］ボルツマン分布ではこのエネルギートラップから現実的な時間内のシミュレーションで抜け出すことが著しく困難である。そのため、多大な量の計算資源と計算時間が消費される。
［阻害要因３］レプリカ交換法などの拡張アンサンブル法を使う事でサンプリング効率を向上させることが可能であるが、温度空間をランダムウォークさせることが難しい。

以上のような阻害要因が存在するため、普通に式（１）のシミュレーションを実行しても最適解にたどり着くことは難しい。そこで、これらの困難を超える方法が提案されており、代表的なものとしてレプリカ交換法がある。レプリカ交換法を採用した場合、最適化装置１００は、異なる温度のシミュレーションを実行しているレプリカ同士を交換しながら全体を最適化することとなる。そのため、計算量は用意するレプリカ数に応じて増大する。また、効率的なサンプリングのためにはすべてのレプリカが低温領域から高温領域までくまなく動き回ることが重要であり、系の自由度が大きくなるとレプリカ数が加速度的に増える。そのうえ、系に相転移が起こると効率的なレプリカ交換が相転移点を境にして阻害され、低温領域と高温領域に二分されてしまう傾向がある。レプリカ数を増やせば相転移が起こり難くなるが、計算量が大きくなり、最適化装置１００として並列計算機を用いることとなる。

そこで、第２の実施の形態では、最適化装置１００は、レプリカに適切な温度設定を与えることで、少ないレプリカ数でも相転移の発生を抑制し、レプリカが温度空間をランダムウォークできるようにする。

次に、レプリカ交換法について詳細に説明する。まず、求解対象のエネルギー関数の構造について簡単に説明する。実用的な問題は多次元であり、求めるエネルギー関数は多数の極小値を持つ。状態変数に応じたエネルギー関数の値の変化を表す図をエネルギーランドスケープという。

図３は、エネルギーランドスケープの一例を示す図である。図３において、横軸は状態変数の値であり、縦軸はエネルギーの値である。多次元空間におけるエネルギー最小値問題の難しさはエネルギー関数の大局的な最小構造にある。たとえば、図３においては大きく見ると３つの大局的なエネルギー最小構造３１〜３３がある。それぞれのエネルギー最小構造３１〜３３内および周囲には、小さいエネルギー障壁（隣接する極小値間を遷移する際のエネルギーの上昇部分）が無数にある。レプリカの設定温度が不適切な場合、初期値によってエネルギー最小構造３１，３３内部でエネルギー最小値を求めてしまい、本当に知りたいエネルギー最小構造３２内のエネルギー最小値を求めることができないことがある。これはエネルギー最小構造３１〜３３間のエネルギー障壁が大きいためである。

レプリカ交換法はこの問題を解決するために開発された方法である。
図４は、レプリカ交換法の概略を示す図である。最適化装置１００は、レプリカ交換法では高温から低温領域までＭ個（Ｍは２以上の整数）の温度列｛Ｔ₀，Ｔ₁，・・・，Ｔ_M｝を用意する。最適化装置１００は、Ｍ個の求解対象（状態変数の組）を生成し、それぞれ温度列内の別の温度により、時間ステップを進行させながらシミュレーションを独立に行う。そして最適化装置１００は、求解対象間で温度を詳細つり合いの原理を満たすように交換する。このようにすると全体として詳細つり合いの原理を満たしつつ、各シミュレーションボックス（求解対象ごとのシミュレーション）は温度空間をランダムウォークすることができる。そして、高温領域を経て高いエネルギーバリアを超える事ができる。

各シミュレーションボックスはレプリカと呼ばれ、２つのレプリカに設定されている温度を一定の基準で交換するところから、このようなシミュレーション方法はレプリカ交換法と呼ばれている。レプリカ間の温度の交換を行うことは、レプリカ交換と呼ばれる。

レプリカ交換法は磁性体分野やタンパク質のフォールディング問題などで成果を上げている方法である。半面、欠点も存在する。以下、レプリカ交換法の欠点について簡単に説明する。

今、ｉ番目とｊ番目との２つのレプリカを考える。このとき逆温度と座標（状態変数の値）の組み合わせを（β_i，｛ｘ｝_i）、（β_j，｛ｘ｝_j）とする。また、それぞれのレプリカの粒子が従う確率分布をπ（β_i，｛ｘ｝_i）、π（β_j，｛ｘ｝_j）とする。この２状態系においてレプリカ交換が起こる前後の状態をＡ、Ｂとして以下のように定義する。
［状態Ａ］ｉ番目のレプリカが（｛Ｘ_i｝，β_i）かつ、ｊ番目のレプリカが（｛Ｘ_j｝，β_j）
［状態Ｂ］ｉ番目のレプリカが（｛Ｘ_j｝，β_i）かつ、ｊ番目のレプリカが（｛Ｘ_i｝，β_j）
熱平衡状態において状態Ａと状態Ｂそれぞれの確率分布Ｐ_A（ｔ）、Ｐ_B（ｔ）は、独立事象の確率の積で表される。すなわち確率分布Ｐ_A（ｔ）、Ｐ_B（ｔ）は、以下の式（５）、式（６）で表される。

従って、状態Ａから状態Ｂへの遷移確率Ｐ_A→Bは以下の式（７）のように記述できる。

ｍｉｎ｛｝は、括弧内の要素の最小値を採ることを示す演算子である。式（７）にボルツマン分布の式を代入すると、以下のようになる。

式（１０）のＥ_iとＥ_jは、それぞれｉ番目のレプリカとｊ番目のレプリカのエネルギーである。各レプリカのエネルギーは、そのときの状態変数を式（１）に代入することで得られる。

レプリカ交換法においては、レプリカ同士の交換についての詳細つり合いが成立し、かつ各レプリカについての詳細つり合いも成立している。なおレプリカ交換は１ステップごとに実施してもよいし、任意のステップごとにレプリカ交換を実施してもよい。ただし、レプリカ交換の頻度によってサンプリング効率が変わる。

このようなレプリカ交換法にも欠点がある。１つ目の欠点は計算量の多さである。計算の実行はＭ個のレプリカを用意した場合、計算量はＭ倍になってしまう。さらに、問題の自由度が大きくなると、計算に必要なレプリカ数も増大してしまう。そのため、通常は並列計算をさせて計算時間を圧縮するが、その代わり要求される計算機の数はＭ倍になる。

２つ目の欠点は温度空間をランダムウォークすることが難しいことである。
図５は、レプリカ交換法におけるレプリカの温度遷移の一例を示す図である。図５には、レプリカ番号＝１〜８までの８つのレプリカのそれぞれにおける、時間ステップの進行に応じた温度（Ｔ）の遷移が示されている。横軸はステップ数を表し、縦軸はＴ（正の実数）である。図５に示されているように、各レプリカは温度空間を等確率でランダムウォークしているとは言えないことが分かる。たとえばレプリカ番号＝１のレプリカは温度が１０以下の領域には殆ど遷移しない。さらにレプリカ番号＝５のレプリカは、ステップ数が２００付近から、Ｔ＝１．０に固定されてしまっている。このように、レプリカ交換法においては、温度空間をランダムウォークさせるのが難しい。

図６は、レプリカ交換法におけるエネルギーの計算結果の一例を示す図である。図６には、レプリカ番号＝１〜８までの８つのレプリカのそれぞれにおける、エネルギー（Ｅ）の計算結果が示されている。横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。

図６に示されているように、殆どのレプリカのエネルギーは極小値に拘束され、脱出できないことが分かる。すなわちランダムウォークができていないと、レプリカのエネルギーは、最小値を含むよりエネルギーの低い値に到達できない。

以下、レプリカ交換の温度最適化方法について詳細に説明する。
まず、レプリカ交換法におけるレプリカの温度を最適化することで最良のサンプリング効率が得られる理由について説明する。レプリカ交換法では、レプリカの温度を一定間隔で設定しても、最良のサンプリング効率は得られない。その理由は、交換確率が、エネルギー差と逆温度差の積の指数関数減衰になるからである。そのため、レプリカの温度が一定間隔の場合、レプリカが高温層と低温層に分離してしまい、高温領域でエネルギーバリアを超えられず、図６に示したようにサンプリング空間が広がらないことが起こる。このように、レプリカの温度が不適切な場合、適切なランダムウォークが起きない。結果としてレプリカ交換法を用いたシミュレーションにおけるサンプリング空間が広がらないことになる。この問題は、求解対象の状態変数の数が少ないうちはあまり問題にならないが、自由度が大きくなると深刻な問題になる。

次に、図７を参照して、レプリカ間の温度差とレプリカ交換確率との関係について説明する。
図７は、温度の異なる２つのレプリカのエネルギーの時間変化の例を示す図である。図７の例では、温度以外は同じ初期条件のレプリカそれぞれについて、エネルギー関数の計算シミュレーションを繰り返し実行したものである。各レプリカは、計算を繰り返すと、熱平衡状態となり、ある一定の範囲内でエネルギー値が変動する。図７において、横軸はステップ数を表し、縦軸はエネルギー（Ｅ）を表す。

図７の左図は、Ｔ＝１００のレプリカ（下の折れ線）とＴ＝３００のレプリカ（上の折れ線）とのエネルギーの時間変化である。この例ではレプリカ間の温度差が大きいためエネルギーの重なり合いが殆どない。そのため、レプリカ交換自体が起こらない。

図７の右図は、Ｔ＝５０のレプリカ（下の折れ線）とＴ＝１００のレプリカ（上の折れ線）のエネルギーの時間変化である。この例ではレプリカ間の温度差が小さいためエネルギーの重なり合いが起こる。エネルギーが重なり合うとレプリカ交換が発生しやすくなる。

このように、レプリカ間の温度差を小さくすれば、レプリカ交換が発生しやすくなる。すべてのレプリカを温度順に並べたときに隣接するレプリカ間の温度差がすべて小さければ、各レプリカは温度空間全体を適切にランダムウォークすることが容易となる。ただし、温度空間全体にわたってレプリカ間の温度差を小さくしようとすると、レプリカの数を増加させることとなる。すると、レプリカ交換法の１つ目の欠点で説明した処理量の増大が顕著となる。

そこで最適化装置１００は、レプリカの温度の適性化を図ることで、レプリカ数の増加を抑止しつつ、レプリカ交換の発生をしやすくし、各レプリカのランダムウォークを実現する。すなわち最適化装置１００は、温度が隣接するレプリカ間でエネルギーの適切な重なりが発生するような温度差となるように、各レプリカの温度を決定する。

たとえば図７の左図に示すように、Ｔ＝１００のレプリカとＴ＝３００のレプリカとは、時間変化に伴うエネルギーの振れ幅が異なる。エネルギーの振れ幅が大きい温度帯域であれば、レプリカ間の温度差を大きくしてもエネルギーの重なり合いが発生し、レプリカ交換が起こりやすい。他方、エネルギーの振れ幅が小さい温度帯域では、レプリカ間の温度差が大きいとエネルギーの重なり合いが発生しない。そこで最適化装置１００は、温度帯域ごとに、その温度帯域におけるレプリカ間の温度差を適切に設定する。これにより、レプリカ数の増加を抑止しながら、温度空間全体にわたってレプリカ交換を起こりやすくする。

次に、最良のサンプリング効率を与えるための温度最適化方法について詳細に説明する。以下の方法において比熱は連続関数であると仮定する。従って、比熱が不連続になるような相転移が起きる場合は対象としない。これはシミュレーションを実施する系において粒子数は常に有限であり、相転移による特異性が現れるのは粒子数無限大の極限であるからである。

なお、レプリカ交換法において最良のサンプリング効率が実現される場合とは、すべてのレプリカが温度空間を上から下までランダムウォークする場合である。このようなランダムウォークを、適切なランダムウォークと呼ぶこととする。

適切なランダムウォークを実現するには、温度順でレプリカを並べたときの隣接するレプリカ間のレプリカ交換の発生確率が、すべて等確率であればよい。従って式（８）においてすべてのレプリカが以下の式（１１）式を満たせばよいことになる。

Δ₀は、レプリカ交換確率係数である。ここでｊ＝ｉ＋１とすると、Ｔ_j＝Ｔ_i+1＝Ｔ_i＋ΔＴ_iとして、式（１２）が得られる。

また、ΔＥは式（１３）のように表記できる。

この場合、レプリカ交換確率係数Δ₀は、以下のように変形される。

ここで、エネルギー差が温度差に対して激しく変化しない場合、あるいは、多数のレプリカがある場合を考えると、次の式（１５）が近似的に成立する。

式（１５）を式（１４）に適用すると、式（１４）は次の式（１６）となる。

式（１６）をΔＴ_iについて解くと、次の式（１７）が得られる。

ここで、Ｃ_Vは定積比熱と呼ばれる量である。以下の説明では、定積比熱Ｃ_Vを単に比熱と呼ぶこともある。またΔＴ_iは、次の式（１８）で表すこともできる。

ある与えられた温度Ｔ_iからレプリカ交換確率を一定にする次の温度Ｔ_i+1は次の式（１９）および式（２０）で与えられる。

つまり、ｉ＋１番目のレプリカの温度Ｔ_i+1は、基準となる温度に対して常に比例係数をかけた形で定義される。次に比熱の振る舞いについて境界条件を与える。比熱は次の境界条件を満たす。

式（２１）は熱力学第三法則からくる要請である。これより温度刻み係数γ（Ｔ）の振る舞いは次の境界条件を満たす。

具体的な問題について比熱Ｃ_V（Ｔ）を計算すると、たとえば図８のようになる。
図８は、比熱の温度変化の一例を示す図である。図８の上段には、温度が０〜１０００００Ｋの範囲の比熱の温度変化が示されている。図８の中段には、温度が０〜２００００Ｋの範囲の比熱の温度変化が示されている。図８の下段には、温度が０〜２００Ｋの範囲の比熱の温度変化が示されている。図８に示す各グラフは、横軸が温度、縦軸が比熱である。

図８の比熱の温度変化によれば、式（２１）、式（２２）に示されている比熱の境界条件が満たされていることが分かる。また、比熱は極小値を２つ、極大値を２つ持つ構造であることが分かる。

図８を参照すると、図５に示したように、レプリカ交換時にレプリカが高温相と低温相に分かれる理由が判明する。図８に示すように、比熱が極大値を持つため、温度Ｔが絶対温度で２０Ｋ＜Ｔ＜４０Ｋの温度帯域では温度増加に伴い比熱が増大しており、増大傾向も急激である。

この比熱が急激に増大する温度帯域では、レプリカ間の温度差によって生じるそれらのエネルギー差も急激に増大する。このため比熱が急激に増大する温度帯域でレプリカに設定する温度を等間隔に取ってしまうと、エネルギー差が過大となった温度帯域において、レプリカ交換確率に影響する式（１１）のレプリカ交換確率係数Δ₀が大きな値になってしまう。そのため、比熱が増大し、エネルギー差が過大となった温度帯域ではレプリカ交換を行う試行は殆ど棄却されてしまう。一方で比熱が急激に増大する領域を外れ、比熱が小さくエネルギー差も少ない温度帯域ではレプリカ交換が行われる。

その結果、比熱が極大点となる温度を境にして実質的にレプリカ交換の起こる温度帯域が分断されてしまい、高温相と低温相に分かれてしまう。そして、温度が上がる方向へのレプリカ交換も、下がる方向へのレプリカ交換も等しい確率で起こらないといけないことも分かる。なぜならば、一方の方向（温度が下がる方向または上がる方向）のみにレプリカ交換が起こりやすくなると確率流が定常状態ではなくなるからである。つまり比熱の極大点がレプリカ交換の確率流の反射壁の役割をしてしまう。

この問題を解決するために最適化装置１００は、たとえば図８の温度で２０Ｋ＜Ｔ＜４０Ｋのように、比熱が急激に増加し、比熱が大きくなった温度帯域では、レプリカに設定する温度刻み幅を小さくする。また最適化装置１００は、比熱が小さい温度帯域ではレプリカに設定する温度刻み幅を大きくする。これにより最適化装置１００は、温度が下がる方向にも上がる方向にも等しい確率でレプリカ交換を起こすようにし、確率の流れを一定に保つことができる。

具体的には最適化装置１００は、設定温度の順にレプリカを並べたときに、隣接するレプリカ間についてレプリカ交換確率係数Δ₀が一定であるという条件を満たすように、レプリカ交換における各レプリカの温度を設定する。この条件を満たすように温度間隔を決定する式が式（２０）である。

図９は、レプリカの温度刻み係数と温度との関係を示す図である。図９では、横軸が温度であり、縦軸が温度刻み係数γ（Ｔ）である。そして温度刻み係数γ（Ｔ）をレプリカ交換確率に関連したパラメータであるレプリカ交換確率係数Δ₀が「０．２」の場合（黒の矩形を結ぶ折れ線）と「０．５」の場合（白丸を結ぶ折れ線）とについて図示している。

図９によれば、比熱が急激に変化する領域では温度刻み係数が小さくなり、温度刻み幅は小さくするのが適切であることが分かる。一方で高温領域では、温度刻み係数がほぼ一定と見なすことができる。このため、高温領域ではレプリカ温度を等間隔に設定してもよいことが分かる。

さらに、温度刻み幅のレプリカ交換確率係数Δ₀の値への依存性も無視できない。なおレプリカ交換確率係数Δ₀は、レプリカ交換法によるシミュレーションを指示するユーザが、シミュレーションの計算量や精度に応じて最適化装置１００に設定する値である。

レプリカ交換確率係数Δ₀が大きくなると、式（８）と式（１１）からレプリカ交換確率は小さくなる。レプリカ交換確率が小さくてよいなら温度刻み幅を大きく取ることができる。この場合、レプリカ交換確率は下がるがレプリカ数は少なくなるため計算量は下がる方向に向かう。他方、レプリカ交換確率を大きくしたいならば、ユーザはレプリカ交換確率係数Δ₀を小さくすればよい。この場合、温度刻み幅は細かくなるが、レプリカ交換確率は大きくなる。そしてレプリカ数は増加し、計算量は増大する方向に向かう。

最大のサンプリング効率の最適化を実現するためには、レプリカ交換確率を増大させ、レプリカ数を増やさずに、交換確率とレプリカ数の両者のトレードオフから決まる最適値を求めるのが適切である。しかし、この最適値は問題依存であるため、問題に応じてユーザが適切に決定することとなる。

次に、最適化装置１００がレプリカ交換法による解求解を行うために有する機能について説明する。
図１０は、最適化装置の機能を示すブロック図である。最適化装置１００は、記憶部１１０と処理部１２０とを有する。記憶部１１０は、たとえばメモリ１０２またはストレージ装置１０３の記憶領域の一部によって実現される。処理部１２０は、たとえばプロセッサ１０１に所定のプログラムを実行させることにより実現される。

記憶部１１０は、温度最適化情報、エネルギー情報、スピン情報、レプリカ情報、温度情報、問題設定情報、ハミルトニアン情報を記憶する。
温度最適化情報は、レプリカの設定温度の最適化に使用する情報である。たとえば温度最適化情報には、レプリカ交換確率係数Δ₀、予備計算ステップ数Ｎ、比熱計算用温度列、最適化する温度列の最小値Ｔ₁’などを含む。レプリカ交換確率係数Δ₀は、レプリカ交換の確率を決定する係数であり、式（１１）を満たす定数である。予備計算ステップ数Ｎは、比熱を計算する温度（比熱計算用の標本温度）でのレプリカのエネルギー計算を計算する際の、シミュレーションの時間進行のステップ数である。比熱計算用温度列は、Ｍ個（Ｍは１以上の整数）の比熱計算用の温度値の配列である。比熱計算用温度列は、第１の実施の形態における複数の標本温度値の一例である。最適化する温度列の最小値Ｔ₁’は、レプリカに設定する温度の最小値である。

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

処理部１２０は、制御部１２１、設定読込部１２２、問題設定部１２３、比熱計算部１２４、温度最適化部１２５、およびエネルギー最適化部１２６を有する。
制御部１２１は、処理部１２０の各部を制御して最小値求解処理を行う。

設定読込部１２２は、記憶部１１０から上記の各種情報を、制御部１２１が理解可能な形式で読み込む。
問題設定部１２３は、状態変数の初期値や、エネルギーの最小値の初期値を計算する。

比熱計算部１２４は、求解問題に対して、複数の温度での比熱を計算する。
温度最適化部１２５は、温度に応じた比熱に基づいて、レプリカに設定する温度の最適化を行う。

エネルギー最適化部１２６は、最適化された温度をレプリカに設定し、レプリカ交換法によるシミュレーションを実行し、エネルギーの最小値を求解する。
なお、図１０に示した各要素間を接続する線は通信経路の一部を示すものであり、図示した通信経路以外の通信経路も設定可能である。また、図１０に示した各要素の機能は、たとえば、その要素に対応するプログラムモジュールをコンピュータに実行させることで実現することができる。

次に、比熱計算部１２４による比熱の計算法について詳細に説明する。比熱は温度の関数として計算することができる。そのため、比熱計算部１２４は、離散的な温度列（温度Ｔ_j（ｊ＝１，２，・・・，Ｍ））に基づいて、式（１）をメトロポリス−ヘイスティングス法によって解く。さらに比熱計算部１２４は、得られたエネルギー列からエネルギーの平均値＜Ｅ_j＞と二乗平均値＜Ｅ²＞とを計算する。そして比熱計算部１２４は、分布がボルツマン分布であるため、次の式（２５）で比熱を計算する。

このように比熱計算部１２４が各温度Ｔ_jにおいて予備計算としてのシミュレーションを短時間行い、比熱を計算することで、図８に示したような温度ごとの比熱が得られる。
以上が比熱の計算方法である。次に温度最適化部１２５による温度の最適化方法について説明する。

温度最適化部１２５は、数値的に比熱が求められると、式（１７）を用いて温度刻み幅ΔＴ_iを計算する。このとき式（１７）は予備計算の時に指定した温度列（温度Ｔ_j（ｊ＝１，２，・・・，Ｍ））以外の温度でも計算できるようにするため、温度最適化部１２５は、補間法などを用いて、比熱を温度に応じた連続変数に拡張する。簡単な例では、温度最適化部１２５は直線補間を用いることができる。たとえば温度最適化部１２５は、温度Ｔ₀がＴ_j＜Ｔ₀＜Ｔ_j+1を満たす場合、たとえば直線補間をする以下の式（２６）で、温度刻み係数γ（Ｔ）を計算する。

式（２６）は直線補間の例であるが、微分まで含め滑らかに連続に補間するような関数を用いて高次の補間法を用いてもよい。
次にエネルギー最小値求解処理の手順について詳細に説明する。

図１１は、エネルギー最小値求解処理の手順の一例を示すフローチャートである。以下、図１１に示す処理をステップ番号に沿って説明する。
［ステップＳ１０１］設定読込部１２２は、記憶部１１０からレプリカ交換確率係数Δ₀を読み込む。レプリカ交換確率係数Δ₀は、ユーザにより予め入力された値である。

レプリカ交換確率係数Δ₀の値が大きいほど温度刻みの幅は大きくなり、レプリカ数は少なくなる。逆にレプリカ交換確率係数Δ₀の値が小さいほど温度刻みの幅は小さくなり、レプリカ数は多くなる。レプリカ交換確率係数Δ₀の値が大きすぎると温度刻みが荒くなりすぎるため、ユーザは、適切なレプリカ数となるように、適切な値のレプリカ交換確率係数Δ₀を設定する。たとえばユーザは、Δ₀＝１．０とする。設定読込部１２２は、読み込んだレプリカ交換確率係数Δ₀を制御部１２１に送信する。

［ステップＳ１０２］設定読込部１２２は、記憶部１１０から予備計算ステップ数Ｎを読み込む。予備計算ステップ数Ｎは、ユーザにより予め入力された値である。
予備計算ステップ数Ｎは大きく取ればとるほどエネルギーの期待値の精度が上がる。他方、大きく取りすぎると計算時間が長くなりすぎる。そのため予備計算ステップ数Ｎは、温度最適化を行うことができる程度に十分小さな値とするのが適切である。たとえばユーザは、予備計算ステップ数Ｎ＝１０⁶回と入力する。この程度の予備計算ステップ数Ｎがあれば、エネルギー関数のエネルギー値が熱平衡状態となるために、熱平衡状態でのエネルギー値を求めることができる。設定読込部１２２は、読み込んだ予備計算ステップ数Ｎを制御部１２１に送信する。

［ステップＳ１０３］設定読込部１２２は、記憶部１１０から、比熱計算用の複数の温度Ｔ_j（ｊ＝１，２，・・・，Ｍ）のリストである比熱計算用温度列を読み込む。比熱計算用の温度Ｔ_jは、ユーザにより予め入力された値である。たとえばユーザは、比熱を計算する温度の離散点を、リストの形式で比熱計算用の温度Ｔ_jとして入力する。設定読込部１２２は、読み込んだ比熱計算の温度Ｔ_jの温度列を制御部１２１に送信する。この際、制御部１２１は、比熱計算用の温度Ｔ_jの個数Ｍを計算する。

［ステップＳ１０４］設定読込部１２２は、記憶部１１０から、最適化後の温度列の個数Ｍ_optを読み込む。最適化後の温度列の個数Ｍ_optは、ユーザが予め入力する値である。設定読込部１２２は、読み込んだ個数Ｍ_optを制御部１２１に送信する。

［ステップＳ１０５］設定読込部１２２は、記憶部１１０から、最適化する温度列の最小値Ｔ₁’の値を読み込む。最適化する温度列の最小値Ｔ₁’は、ユーザが予め入力する値である。設定読込部１２２は、読み込んだ最小値Ｔ₁’の値を制御部１２１に送信する。

［ステップＳ１０６］制御部１２１は、比熱計算用温度列に示される温度における比熱を、比熱計算部１２４に計算させる。たとえば制御部１２１は、比熱計算部１２４へ、レプリカ交換確率係数Δ₀、予備計算ステップ数Ｎ、比熱計算用温度列（温度Ｔ_j（ｊ＝１，２，・・・，Ｍ））、最適化する温度列の最小値Ｔ₁’などの情報を送信する。比熱計算部１２４は、記憶部１１０からレプリカのエネルギー計算に用いる情報を読み込み、比熱計算処理を行う。比熱計算処理の詳細は後述する（図１２参照）。

［ステップＳ１０７］制御部１２１は、レプリカに設定する温度の最適化処理を、温度最適化部１２５に実行させる。たとえば制御部１２１は、比熱計算部１２４が計算した、比熱計算用の温度ごとの比熱を、温度最適化部１２５に送信する。温度最適化部１２５は、比熱計算用の温度ごとの比熱に基づいて、レプリカに設定する最適な温度を計算する。温度最適化処理の詳細は後述する（図１３参照）。

［ステップＳ１０８］制御部１２１は、エネルギー最適化部１２６にエネルギー最適化処理を実行させる。たとえば制御部１２１は、温度最適化部１２５が計算した温度の値を、エネルギー最適化部１２６に送信する。エネルギー最適化部１２６は、取得した温度をレプリカの温度に設定し、レプリカ交換法による求解問題のシミュレーションを実行し、エネルギーの最小値を算出する。

次に、比熱計算処理について詳細に説明する。
図１２は、比熱計算処理の手順の一例を示すフローチャートである。以下、図１２に示す処理をステップ番号に沿って説明する。

［ステップＳ１２１］比熱計算部１２４は、比熱計算用温度列における各温度の順番を示す変数ｊの値を１からＭまで１ずつ増加させていき、変数ｊに設定された値それぞれについて、ステップＳ１２２〜Ｓ１２３の処理を実行する。

［ステップＳ１２２］比熱計算部１２４は、エネルギー最適化部１２６に、温度Ｔ_jにおけるイジングモデルのサンプリング計算をＮ回実行させる。たとえば比熱計算部１２４は、制御部１２１に比熱計算用の温度Ｔ_jでのＮ回のサンプリング計算を依頼する。すると制御部１２１は、エネルギー最適化部１２６に対して、サンプリング計算を指示する。

エネルギー最適化部１２６は、記憶部１１０から、スピン情報、レプリカ情報、温度情報、問題設定情報、ハミルトニアン情報などを取得する。そしてエネルギー最適化部１２６は、式（１）に従って、メトロポリス−ヘイスティングス法によってイジングモデルのエネルギーの計算をＮ回行う。このようなサンプリング計算によるエネルギー値は、各比熱計算用の温度Ｔ_jで決まる熱平衡状態に緩和する。エネルギー最適化部１２６は、エネルギーの計算を行うごとに、算出したエネルギー値を記憶部１１０に格納する。なお、格納されたエネルギー値の集合は、図１に示す熱平衡状態情報１１ｂの一例である。エネルギー最適化部１２６は、Ｎ回の計算が終了すると、サンプリング計算の終了を制御部１２１に通知する。制御部１２１は、比熱計算部１２４にサンプリング計算の完了を通知する。

［ステップＳ１２３］比熱計算部１２４は、平衡状態になったエネルギーの平均値＜Ｅ_j＞および二乗平均値＜Ｅ_j ²＞を計算し、式（２５）により比熱Ｃ_v（Ｔ_j）を計算する。そして比熱計算部１２４は、式（２０）により温度刻み係数γ（Ｔ_j）の値を計算する。

［ステップＳ１２４］比熱計算部１２４は、１からＭまでの変数ｊの値それぞれについてのステップＳ１２２〜Ｓ１２３の処理が終了すると、比熱計算処理を終了する。その後、比熱計算部１２４は、比熱計算用の温度Ｔ_jごとの比熱Ｃ_v（Ｔ_j）と温度刻み係数γ（Ｔ_j）とを、計算結果として制御部１２１に送信する。

比熱計算処理が終了すると、制御部１２１は、比熱計算用の温度Ｔ_jごとの比熱Ｃ_v（Ｔ_j）と温度刻み係数γ（Ｔ_j）とを温度最適化部１２５に送信し、温度最適化部１２５に温度最適化処理を実行させる。

図１３は、温度最適化処理の手順の一例を示すフローチャートである。以下、図１３に示す処理をステップ番号に沿って説明する。
［ステップＳ１３１］温度最適化部１２５は、最適化後の温度の値の順番を示す変数ｉの値を１からＭ_optまで１ずつ増加させていき、変数ｉに設定された値それぞれについて、ステップＳ１３２〜Ｓ１３３の処理を実行する。

［ステップＳ１３２］温度最適化部１２５は、Ｍ個の（Ｔ_j，γ（Ｔ_j））の組からγ（Ｔ_i）の値を算出する。たとえば温度最適化部１２５は、Ｍ個の（Ｔ_j，γ（Ｔ_j））の組に基づいて式（２６）を用いて補間計算をして、γ（Ｔ_i）の値を算出する。

［ステップＳ１３３］温度最適化部１２５は、次の最適化された温度Ｔ_i+1’を算出する。次の最適化された温度Ｔ_i+1’は、式（１９）に基づいて、Ｔ_i+1’＝γ（Ｔ_i’）Ｔ_i’と計算できる。

［ステップＳ１３４］温度最適化部１２５は、１からＭ_optまでの変数ｉの値それぞれについてのステップＳ１３２〜Ｓ１３３の処理が終了すると、温度最適化処理を終了する。その後、温度最適化部１２５は、最適化された温度Ｔ_i’（ｉ＝１，２，・・・，Ｍ_opt）を、計算結果として制御部１２１に送信する。

このような温度最適化処理により、低い方の温度から順に最適化された温度が決定される。その後、エネルギー最適化部１２６により、最適化した温度を各レプリカの温度に設定して、レプリカ交換法によるシミュレーションが実行される。最適化された温度は、レプリカ交換確率が一定になるような温度刻み幅となっている。そのためエネルギー最適化部１２６によるシミュレーションでは、各レプリカを、温度空間内を適切にランダムウォークさせることができる。その結果、ヒューリスティックに低いエネルギー最適値を得ることができる。

以下、温度の最適化の具体例について説明する。
例題として１６都市からなる巡回セールスマン問題を想定する。巡回セールスマン問題は、都市の集合と都市間の移動コストが与えられたとき、すべての都市を１回だけ巡回するための、総移動コストが最小となる経路を求める組合せ最適化問題である。最適化装置１００が巡回セールスマン問題において適切な解を求解することができれば、最適化装置１００を用いて、配送ルート最適化などの社会科学上の様々な問題に対する適切な解が得られる。

最適化装置１００の比熱計算部１２４は、式（２０）を巡回セールスマン問題に適用することで、図８に示すような比熱を算出する。比熱計算部１２４は、算出した比熱から式（２０）を用いて温度刻み係数γ（Ｔ）を算出する。たとえば、最低温度をＴ₀＝１．０とし、レプリカ交換確率を「０．９」とするようなレプリカ交換確率係数をΔ₀≒０．０４５８とする。この条件で上記の巡回セールスマン問題の温度刻み係数γ（Ｔ）を計算すると、図９のΔ₀＝０．５の場合とほぼ同じ温度刻み係数γ（Ｔ）が得られる。

ここで、レプリカ数を「５」（Ｍ_opt＝５）として最適化した場合、以下のような温度列が得られる。
＜最適化した温度列（最適化例）＞
｛Ｔ₀，Ｔ₁，Ｔ₂，Ｔ₃，Ｔ₄｝＝｛１．０，５．０５５，１３．３４９，２８．３５４，３２．６７７｝
ここで比較例として、温度最適化を行わない場合にレプリカに設定する温度列を以下に示す。
＜最適化しない温度列（第１の比較例）＞
｛Ｔ₀，Ｔ₁，Ｔ₂，Ｔ₃，Ｔ₄｝＝｛１．０，２０．０，４０．０，６０．０，８０．０｝
比較例は、ユーザが比熱に関する予備知識を持たず、温度の選択をする基準をもたない場合を想定している。ユーザに予備知識がなければ、ユーザは、比較例に示すような一定間隔の温度列を設定されるものと推定できる。

また、温度最適化によるエネルギー最小値候補を得られる効率が、他のサンプリング手法を用いた場合に比べて有意に改善するかどうかを確認するため、ダイナミックオフセット法と比較する。ダイナミックオフセット法とは、多変数関数の最小値問題において、解（状態変数の組）がエネルギーの極小値にはまり脱出できなくなった場合に使う方法である。具体的には、ダイナミックオフセット法は、エネルギーの値に一定の値を加算していくことで、解を強制的にエネルギー極小値から脱出させるのである。ダイナミックオフセット法を使うことで、解がエネルギー極小値にトラップされることがなくなり、エネルギーの最小値候補を効率的に得られることが知られている。

図１４は、計算ステップ数を増やした場合と温度最適化を行った場合との解を比較した例を示す図である。第１のケース４１は、上記の第１の比較例に示される温度列を使用し、計算回数は２００万回として最小値求解を行った例である。第２のケース４２は、上記の第１の比較例に示される温度列を使用し、計算回数は８００万回として最小値求解を行った例である。第３のケース４３は、上記の最適化例に示される温度列を使用し、計算回数は２００万回として最小値求解を行った例である。

図１４の横軸は、各ケースで得られた解についてエネルギーが小さい順に順位付けしたときの、同一ケース内での順位（ランキング）である。縦軸は、各ランキングの解のエネルギーの値である。なお何れのケースでもダイナミックオフセット法を併用している。これは極小値にはまらないような状況下で、温度最適化の効果のみを定量的に示すためである。

図１４に示されるように、温度最適化を行った第３のケース４３では、同じ計算回数の温度最適化を行っていない第１のケース４１よりもエネルギーの低い解が得られている。温度最適化を行わずに計算回数を４倍にした第２のケース４２では、第１のケースよりもエネルギーの低い解が得られているが、第３のケース４３程ではない。すなわち、温度最適化を行わずに第３のケース４３と同程度の解を得るには、計算回数を４倍以上とすることが求められる。換言すると、温度最適化をすることで計算効率が４倍以上となっている。これらのことから、以下のような結論が得られる。

＜結論１＞
温度最適化することで、シミュレーションで得られるエネルギー最小値と次善の解（良解）のエネルギー値は有意に下がる。

＜結論２＞
温度最適化をすることで、エネルギー最小値と次善の解（良解）を得るための計算量は温度最適化をしない場合よりも少なくなる。

次に、レプリカ数を増やした場合と温度最適化を行った場合とにおける最小値求解の結果を比較する。レプリカ数を増やした場合の例として、以下のような最適化しない温度列を用いるものとする。
＜最適化しない温度列（第２の比較例）＞
｛Ｔ₀，Ｔ₁，Ｔ₂，Ｔ₃，Ｔ₄，Ｔ₅，Ｔ₆，Ｔ₇，Ｔ₈，Ｔ₉，Ｔ₁₀，Ｔ₁₁，Ｔ₁₂，Ｔ₁₃，Ｔ₁₄，Ｔ₁₅，Ｔ₁₆｝＝｛１．０，５．０，１０．０，１５．０，２０．０，２５．０，３０．０，３５．０，４０．０，４５．０，５０．０，５５．０，６０．０，６５．０，７０．０，７５．０，８０．０｝
第２の比較例は、レプリカ数を１７個にした場合の各レプリカに設定する温度列である。温度列内の各温度は、第１の比較例と同様に等間隔の温度となっている。

図１５は、レプリカ数を増やした場合と温度最適化を行った場合との解を比較した例を示す図である。第４のケース５１は、上記の第１の比較例に示される温度列を使用し、５つのレプリカにより最小値求解を行った例である。第５のケース５２は、上記の第２の比較例に示される温度列を使用し、１７個のレプリカにより、計算回数は２００万回として最小値求解を行った例である。第６のケース５３は、上記の最適化例に示される温度列を使用し、５つのレプリカにより、計算回数は２００万回として最小値求解を行った例である。なお、いずれのケースにおいても、計算回数は２００万回であり、レプリカの交換のための試行頻度は、１０回の計算に１回である。

図１５の横軸は、各ケースで得られた解についてエネルギーが小さい順に順位付けしたときの、同一ケース内での順位（ランキング）である。縦軸は、各ランキングの解のエネルギーの値である。なお何れのケースでもダイナミックオフセット法を併用している。

図１５に示されるように、温度最適化を行った第６のケース５３では、レプリカ数の温度最適化を行っていない第４のケース５１よりもエネルギーの低い解が得られている。温度最適化を行わずにレプリカ数を１７に増やした第５のケース５２では、第４のケース５１よりもエネルギーの低い解が得られているが、第６のケース５３程ではない。すなわち、温度最適化を行った第６のケース５３では、５つのレプリカの系でも、レプリカ数を１７個まで増やした第５のケース５２より良いエネルギー最小値候補が得られている。これらのことから、以下のような結論が得られる。

＜結論３＞
レプリカ交換温度を最適化することで、必要レプリカ数を削減することができる。
なお、レプリカの数は計算量にほぼ比例する。そのため、たとえばレプリカ数を１７個から５個に削減することができれば、計算量は３分の１以下に削減できる。

以上より、レプリカ交換の温度の最適化が可能な最適化装置１００を用いることで、以下のような技術的な効果が得られる。
１．レプリカ交換法による最小値求解におけるレプリカの適切なランダムウォークによる精度の高い解の取得。
２．レプリカ交換法による最小値求解に使用する計算資源の削減。
３．レプリカ交換法による最小値求解に要する計算時間の削減。
４．レプリカ交換法による最小値求解に最適なレプリカの温度列の自動決定によるユーザの労力削減。

〔第３の実施の形態〕
次に第３の実施の形態について説明する。第３の実施の形態は、レプリカに設定する温度の精度を向上させるものである。すなわち第３の実施の形態では、最適化装置１００において、比熱計算を繰り返し行うことにより、レプリカに設定する温度の精度を向上させる。

最適化装置１００は、比熱の計算からレプリカ交換法においてランダムウォークを実現するために適切な温度列を自動決定する。一方で最適化された温度の精度は比熱を計算するために実行した予備計算に依存する。

予備計算時には温度は最適化されていないため、予備計算（図１２のステップＳ１２２）に使用したサンプリングの温度は最良ではない。そのため、比熱の精度も近似的な見積もりになっている。比熱の精度を向上させるためには、予備計算において、適切なサンプリング温度を設定することが重要である。予備計算における適切なサンプリング温度は、温度最適化における温度刻み係数γ（Ｔ）の算出（図１３のステップＳ１３２）において比熱を正しく算出可能な温度である。

たとえば温度最適化処理においてレプリカの温度として最適化した温度を用いて予備計算を行うことで、温度刻み係数γ（Ｔ）算出時に比熱を補間で求めたことによる誤差を最小限に抑えることができる。すなわち、最適化された温度列を比熱計算用温度列として予備計算を行うことで、温度刻み係数γ（Ｔ）の計算に使用する比熱の精度が向上し、最適化された温度の精度も向上する。

図１６は、第３の実施の形態におけるエネルギー最小値求解処理の手順の一例を示すフローチャートである。図１６に示す処理のうちステップＳ２０１〜Ｓ２０７，Ｓ２１０は、それぞれ図１１に示した第２の実施の形態の処理のステップＳ１０１〜Ｓ１０８の処理と同様である。そこで以下に、第２の実施の形態と異なるステップＳ２０８〜Ｓ２０９の処理について説明する。

［ステップＳ２０８］制御部１２１は、ステップＳ２０６〜Ｓ２０７によるレプリカの温度の決定処理を所定回数繰り返したか否かを判断する。制御部１２１は、所定回数繰り返していれば処理をステップＳ２１０に進める。また制御部１２１は、所定回数繰り返していなければ、処理をステップＳ２０９に進める。

［ステップＳ２０９］制御部１２１は、比熱計算用温度列に、ステップＳ２０７で最適化した温度列を設定する。その後、制御部１２１は、処理をステップＳ２０６に進める。
このように、最適化装置１００は、最適化された温度列を用いた比熱の再評価と、精度を向上させた比熱による温度の最適化を繰り返し実行することで、温度最適化の精度を向上させることができる。レプリカの温度の最適化の精度が向上すれば、ランダムウォークが行われる可能性が高くなり、エネルギー最小値候補（良解）の精度も向上する。

〔その他の実施の形態〕
最適化装置１００は、温度を最適化したレプリカ交換法にダイナミックオフセット法を組み合わせてシミュレーションを行うことができる。ダイナミックオフセット法を組み合わせることで、さらにヒューリスティックに低いエネルギーを得ることができる。

なお、第２の実施の形態に示したように、レプリカ交換確率はエネルギーと温度列のみで決まる。従って最適化装置１００で実現したレプリカ交換法における温度最適化は、イジングモデルや巡回セールスマン問題以外の様々な最適化問題の求解にも利用することができる。

たとえば、計算化学の分野では分子動力学法が頻繁に使われる。分子動力学法でもレプリカ交換法が利用されるが、レプリカ交換法のレプリカの温度の最適化は、イジングモデルと同様に重要である。レプリカ交換確率は温度とエネルギーで決まり、状態変数が離散的であるか連続的であるかに依存しない。そのためユーザが各温度において、各時刻におけるエネルギー値のサンプリング情報さえ最適化装置１００に与えれば、レプリカ交換を最適化する温度列を、イジングモデルの場合と同じ論理で求めることが可能である。

図１７は、レプリカ交換を最適化する温度列を求めるためのサンプリング情報の一例を示す図である。サンプリング情報６０には、サンプリングのためのシミュレーションのステップ数ごとに、レプリカ温度とエネルギーとが設定されている。このようなサンプリング情報６０があれば、最適化装置１００は、温度ごとの比熱を求め、レプリカ交換法におけるレプリカの温度最適化を行うことができる。なお、サンプリング情報６０は、図１に示した第１の実施の形態における熱平衡状態情報１１ｂの一例である。

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

１０ 最適化装置
１１ 記憶部
１１ａ 求解問題情報
１１ｂ 熱平衡状態情報
１２ 処理部

Claims (10)

1. 温度変化に応じたエネルギーの変化度合いを示す比熱と温度との関係を示す比熱情報を取得し、前記比熱情報に基づいて、レプリカ交換法でエネルギー関数の最適値を求解する際に複数のレプリカそれぞれに設定する複数の温度値を決定する処理部、
   を有する最適化装置。
2. 前記エネルギー関数の最適値を求める求解問題における前記エネルギー関数の複数の標本温度値それぞれでの熱平衡状態のエネルギー値を示す熱平衡状態情報を記憶する記憶部をさらに有し、
   前記処理部は、前記熱平衡状態情報に基づいて、前記比熱情報を算出する、
   請求項１記載の最適化装置。
3. 前記処理部は、前記複数の標本温度値を値の大きさで並べたときに隣接する２つの標本温度値それぞれに対応する前記比熱の間を補完する補間式を、前記比熱情報として生成する、
   請求項２に記載の最適化装置。
4. 前記処理部は、前記複数の標本温度値それぞれにおいて、前記エネルギー関数の計算を繰り返し、熱平衡状態になった後のエネルギー値を前記熱平衡状態情報として前記記憶部に格納する、
   請求項２または３に記載の最適化装置。
5. 前記処理部は、前記複数の温度値を決定後、前記複数の標本温度値を前記温度値と同じ値に変更し、変更後の前記複数の標本温度値それぞれにおいて、前記エネルギー関数の計算を繰り返し、熱平衡状態になった後のエネルギー値で前記熱平衡状態情報を更新する、
   請求項４記載の最適化装置。
6. 前記処理部は、第１の温度帯域における温度値同士の間隔が、前記第１の温度帯域よりも比熱が小さい第２の温度帯域における温度値同士の間隔よりも狭くなるように、前記複数の温度値を決定する、
   請求項１ないし５のいずれかに記載の最適化装置。
7. 前記処理部は、前記比熱が大きいほど値が小さくなる温度刻み幅算出式に基づいて、温度が低い方からｍ番目（ｍは１以上の整数）の温度値における前記比熱に応じた温度刻み幅を算出し、前記ｍ番目の温度値よりも算出した前記温度刻み幅だけ高い温度を、温度が低い方からｍ＋１番目の温度値に決定する、
   請求項６記載の最適化装置。
8. 前記処理部は、さらに、
   前記複数のレプリカそれぞれに前記複数の温度値を設定し、前記レプリカ交換法により前記エネルギー関数の最適値を求解する、
   請求項１ないし７のいずれかに記載の最適化装置。
9. コンピュータが、
   温度変化に応じたエネルギーの変化度合いを示す比熱と温度との関係を示す比熱情報を取得し、
   前記比熱情報に基づいて、レプリカ交換法でエネルギー関数の最適値を求解する際に複数のレプリカそれぞれに設定する複数の温度値を決定する、
   最適化方法。
10. コンピュータに、
    温度変化に応じたエネルギーの変化度合いを示す比熱と温度との関係を示す比熱情報を取得し、
    前記比熱情報に基づいて、レプリカ交換法でエネルギー関数の最適値を求解する際に複数のレプリカそれぞれに設定する複数の温度値を決定する、
    処理を実行させる最適化プログラム。
    

Images