JP2022161128A - プログラム、データ処理方法及びデータ処理装置 - Google Patents

Abstract

【課題】組合せ最適化問題の計算時間を短縮する。【解決手段】記憶部１１は、乱数値に基づく対数計算により得られるノイズ関数の上下限値を、乱数値の範囲を分割した各区間について記憶する。処理部１２は、記憶部１１から、イジング型の評価関数に含まれる各状態変数について生成した乱数値が属する区間の上下限値を各状態変数について取得し、１つの状態変数の値が変化した場合の評価関数の値の変化量と０のうちの大きい方の値と温度パラメータとの積と、その状態変数について取得した下限値との和である下限評価値と、その積とその状態変数について取得した上限値との和である上限評価値とを各状態変数について計算し、上下限評価値に基づいて各状態変数の間の評価値の大小関係を判定することで評価値が最小の状態変数を決定し、決定した状態変数の値を変える、処理を繰り返すことで、評価関数の値が極小または極大になる複数の状態変数の値を探索する。【選択図】図１

Description

本発明は、プログラム、データ処理方法及びデータ処理装置に関する。

ノイマン型コンピュータが不得意とする大規模な離散最適化問題を計算するデータ処理装置として、イジング型の評価関数（エネルギー関数などとも呼ばれる）を用いたイジング装置（ボルツマンマシンとも呼ばれる）がある（たとえば、特許文献１参照）。

イジング装置は、組合せ最適化問題を磁性体のスピンの振る舞いを表すイジングモデルに変換する。そして、イジング装置は、疑似焼き鈍し法やレプリカ交換法などのマルコフ連鎖モンテカルロ法により、イジング型の評価関数の値（イジングモデルのエネルギーに相当する）が極小になるイジングモデルの状態の探索を行う。評価関数の極小値のうちの最小値になる状態が最適解となる。なお、イジング装置は、評価関数の符号を変えれば、評価関数の値が極大になる状態を探索することもできる。以下、マルコフ連鎖モンテカルロ法を、ＭＣＭＣ（Markov-Chain Monte Carlo）法と略す。イジングモデルの状態は、複数の状態変数の値の組合せにより表現できる。各状態変数の値として、０または１を用いることができる。

イジング型の評価関数は、たとえば、以下の式（１）で定義される。

右辺の１項目は、イジングモデルの全状態変数の全組合せについて、漏れと重複なく、２つの状態変数の値（０または１）と重み値（２つの状態変数の間の相互作用の強さを表す）との積を積算したものである。ｘ_ｉは、識別番号がｉの状態変数、ｘ_ｊは、識別番号がｊの状態変数であり、Ｗ_ｉｊは、識別番号がｉとｊの状態変数間の相互作用の大きさを示す重み値である。右辺の２項目は、各識別番号についてのバイアス係数と状態変数との積の総和を求めたものである。ｂ_ｉは、識別番号＝ｉについてのバイアス係数を示している。ｃは定数である。

また、ｘ_ｉの値の変化に伴うエネルギーの変化量（ΔＥ_ｉ）は、以下の式（２）で表される。

式（２）において、ｘ_ｉが１から０に変化するとき、Δｘ_ｉは－１となり、状態変数ｘ_ｉが０から１に変化するとき、Δｘ_ｉは１となる。
イジング装置は、解が局所解に陥ることを防ぐために、エネルギーが上がる状態遷移についても所定の確率で受け入れる。ｘ_ｉの値の変化を受け入れる確率として、たとえば、確率Ａ_ｉ＝ｍｉｎ（１，ｅｘｐ（－ΔＥ_ｉ／Ｔ））を用いることができる。Ｔは温度パラメータである。イジング装置は、たとえば、ｘ_ｉ（ｉ＝１，２，…，Ｎ）をシーケンシャルまたはランダムに選択し、０＜ｕ［ｉ］＜１の一様乱数ｕ［ｉ］を用いて、ｕ［ｉ］＜Ａ_ｉの場合に、ｘ_ｉの値を変化させ、それ以外の場合にはｘ_ｉの値を変化させない。

しかし、上記の確率を用いた場合、状態遷移がいつまでも許容されず、同じ状態が繰り返される可能性がある。
そこで、各試行において状態遷移を発生させる棄却フリー（リジェクション－フリー）ＭＣＭＣ法と呼ばれる手法が用いられる場合がある。リジェクション－フリーＭＣＭＣ法では、Ｎ個の状態変数のうち、以下の式（３）に示す重み付き確率Ｐ_ｉに基づいて、各試行において値を変化させる状態変数が１つ選択される。これにより、通常のＭＣＭＣ法を行った場合と等価な結果を少ない試行回数で得ることができる。

重み付き確率Ｐ_ｉに基づいて値を変化させる状態変数を１つ選択するアルゴリズムとして、以下に示すような重み付けサンプリングのアルゴリズム（たとえば、非特許文献１参照）を用いることができる。

そのアルゴリズムは、ｘ_１～ｘ_ＮのそれぞれについてのＡ_１～Ａ_Ｎに基づいて、ｋ_ｉ＝ｕ［ｉ］^{（１／Ａｉ）}が最大となるｉを識別番号としてもつｘ_ｉを選択するという手法である。ここで、Ａ_ｉ＝ｍｉｎ（１，ｅｘｐ（－ΔＥ_ｉ／Ｔ））であるから、ｋ_ｉ＝ｕ［ｉ］^{（１／Ａｉ）}が最大となるｉという条件は、以下の式（４）で表せる。

式（４）において、ノイズ関数ｌｏｇ（－ｌｏｇ（ｕ［ｉ］））に、温度パラメータ（Ｔ）を乗じた値は熱ノイズに相当する。

特開２０２０－１４０６３１号公報

P.S. Efraimidis, P.G. Spirakis, "Weighted random sampling with a reservoir", Information Processing Letters. 97, 181-185, 2006

上記のように重み付き確率に基づいて値を変化させる状態変数を選択するためには、ノイズ関数を計算することになる。ノイズ関数の計算は対数計算を含むため、計算に時間がかかり、それに伴って組合せ最適化問題の計算時間が長くなる可能性がある。

１つの側面では、本発明は、組合せ最適化問題の計算時間を短縮するプログラム、データ処理方法及びデータ処理装置を提供することを目的とする。

１つの実施態様では、複数の状態変数を含むイジング型の評価関数の値が極小または極大となる前記複数の状態変数の値を、乱数値を用いて探索する処理をコンピュータに実行させるプログラムであって、前記乱数値の範囲を分割した複数の区間のそれぞれについて、前記乱数値に基づく対数計算により得られるノイズ関数の下限値と上限値を記憶した記憶部から、前記評価関数に含まれる前記複数の状態変数のそれぞれについて生成した前記乱数値が属する区間の前記下限値と前記上限値を、前記複数の状態変数のそれぞれについて取得し、前記複数の状態変数のうちの１つの状態変数の値が変化した場合の前記評価関数の値の変化量と０のうちの大きい方の値と温度パラメータとの積と、前記状態変数について取得した前記下限値との和である下限評価値と、前記積と前記状態変数について取得した前記上限値との和である上限評価値とを、前記複数の状態変数のそれぞれについて計算し、前記下限評価値と前記上限評価値に基づいて、前記複数の状態変数のそれぞれの間の、前記積と前記ノイズ関数の和による評価値の大小関係を判定することで、前記評価値が最小となる更新対象の状態変数を決定し、決定した前記更新対象の状態変数の値を変える、処理を繰り返すことで、前記評価関数の値が極小または極大になる前記複数の状態変数の値を探索する、処理をコンピュータに実行させるプログラムが提供される。

また、１つの実施態様では、データ処理方法が提供される。
また、１つの実施態様では、データ処理装置が提供される。

１つの側面では、本発明は、組合せ最適化問題の計算時間を短縮できる。

第１の実施の形態のデータ処理装置及びデータ処理方法の一例を示す図である。 ノイズ関数の一例を示す図である。 第２の実施の形態のデータ処理装置のハードウェア例を示すブロック図である。 データ処理装置の機能例を示すブロック図である。 ノイズ関数上限下限テーブルの一例を示す図である。 第２の実施の形態のデータ処理方法の一例の流れを示すフローチャートである。 探索処理の一例の流れを示すフローチャートである。 ｉ_ｍｉｎの選択処理の一例の流れを示すフローチャートである。 評価値とノイズ関数とβΔＥ_ｉの計算例を示す図である。 第３の実施の形態のデータ処理装置の一例を示す図である。

以下、発明を実施するための形態を、図面を参照しつつ説明する。
前述のリジェクション－フリーＭＣＭＣ法では、各試行において、ｎ個の状態変数のうち、式（４）の関係を満たすｉで示される状態変数の値を変更していくことで、評価関数の値（エネルギー）を最小にするｎ個の状態変数の値（最適解）の探索が行われる。

ところで、上記のｉを決定することは、式（４）のｍａｘ（０，ΔＥ_ｉ）＋Ｔｌｏｇ（－ｌｏｇ（ｕ［ｉ］））において、β＝１／Ｔとして、ｆ（ｉ）＝β×ｍａｘ（０，ΔＥ_ｉ）＋ｌｏｇ（－ｌｏｇ（ｕ［ｉ］））を最小にするｉを決定することに等しい。以下では、Ｔの代わりに、Ｔの逆数であるβを温度パラメータと呼ぶことにする。また、以下、ｆ（ｉ）を評価値と呼ぶ場合もある。

以下に示すデータ処理装置及びデータ処理方法は、このような評価値を最小にするｉで示される状態変数の値を変更していくことで、エネルギーを最小にする最適解の探索を行うものである。

（第１の実施の形態）
図１は、第１の実施の形態のデータ処理装置及びデータ処理方法の一例を示す図である。

データ処理装置１０は、たとえば、コンピュータであり、記憶部１１、処理部１２を有する。
記憶部１１は、たとえば、ＤＲＡＭ（Dynamic Random Access Memory）などの電子回路である揮発性の記憶装置、または、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの電子回路である不揮発性の記憶装置である。

記憶部１１は、ノイズ関数上限下限テーブル１１ａを記憶する。
ノイズ関数は、前述のようにｌｏｇ（－ｌｏｇ（ｕ［ｉ］））と表すことができ、ｕ［ｉ］は、０＜ｕ［ｉ］＜１の一様乱数である。

図２は、ノイズ関数の一例を示す図である。横軸ｘはｕ［ｉ］を表し、縦軸ｙはノイズ関数の値の大きさを表す。
ｕ［ｉ］は疑似乱数ｒ［ｉ］を用いて生成することができる。ｒ［ｉ］が３２ビットの整数の疑似乱数であるものとした場合、ｕ［ｉ］＝ｒ［ｉ］／２^３２である。以下の例では、ｒ［ｉ］を乱数値として説明するが、ｕ［ｉ］を直接、乱数値として扱ってもよい。

図１には、ノイズ関数上限下限テーブル１１ａの一例が示されている。
ノイズ関数上限下限テーブル１１ａは、乱数値に基づく対数計算により得られるノイズ関数の下限値（ｌｂ［ｋ］）と上限値（ｕｂ［ｋ］）を、乱数値の範囲を分割した複数の区間のそれぞれについて含む。

ｋは、各区間の識別番号であり、ｒ［ｉ］がとり得る値の範囲をＫ個に分割した場合、図１のように、ｋ＝１～Ｋである。たとえば、ｋ＝１の区間は、ｒ［ｉ］の下限が１、上限が１×Ｄ－１の区間である。ｒ［ｉ］が３２ビットの整数の疑似乱数とした場合、Ｄ＝２^３２／Ｋである。

ｌｂ［ｋ］は、各区間におけるｒ［ｉ］の上限から求めることができ、ｕｂ［ｋ］は、各区間におけるｒ［ｉ］の下限から求めることができる。たとえば、上記のようにｕ［ｉ］＝ｒ［ｉ］／２^３２で、Ｋ＝２５６の場合、ｌｂ［１］＝ｌｏｇ（－ｌｏｇ（Ｄ／２^３２））＝１．７１２９となり、ｕｂ［１］＝ｌｏｇ（－ｌｏｇ（１／２^３２））＝３．０９９となる。

前述のｆ（ｉ）は、ｒ［ｉ］が区間ｋに属する場合、ｌｂ［ｋ］、ｕｂ［ｋ］を用いて表されるｇ（ｉ）＝β×ｍａｘ（０，ΔＥ_ｉ）＋ｌｂ［ｋ］、ｈ（ｉ）＝β×ｍａｘ（０，ΔＥ_ｉ）＋ｕｂ［ｋ］に対して、ｇ（ｉ）≦ｆ（ｉ）≦ｈ（ｉ）という関係を満たす。以下、ｇ（ｉ）を下限評価値、ｈ（ｉ）を上限評価値という場合もある。

なお、記憶部１１は、処理部１２にデータ処理処理を実行させるプログラムを記憶する。また記憶部１１は、たとえば外部から入力された問題情報（式（１）のＷ_ｉｊ、ｂ_ｉ、ｃ）や計算条件（たとえばβを変化させるための情報、最大試行回数（試行回数の上限値）など）、ｎ個の状態変数の値、エネルギーの計算結果など各種の情報を記憶してもよい。

処理部１２は、たとえば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、ＤＳＰ（Digital Signal Processor）などのハードウェアであるプロセッサにより実現できる。また、処理部１２は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの電子回路により実現されるようにしてもよい。処理部１２は、記憶部１１に記憶されたプログラム（プログラム）を実行して、以下に示すような処理をデータ処理装置１０に行わせる。なお、処理部１２は、複数のプロセッサの集合であってもよい。

処理部１２は、ｆ（ｉ）を最小にするｉで示される状態変数の値を変更していくことで、エネルギーが極小となる解の探索を行う。評価関数の極小値のうちの最小値になる割当状態が最適解となる。なお、式（１）に示した評価関数の符号を変えれば、処理部１２は、評価関数の値が極大になる解を探索することもできる（この場合、最大値が最適解となる）。

ｆ（ｉ）＞ｆ（ｊ＝ｉ＋１）であるか否かは、上記のｆ（ｉ）とｆ（ｊ）についての下限評価値と上限評価値（前述のｌｂ［ｋ］とｕｂ［ｋ］を用いて得られる）とを用いて判定できる場合がある。すなわち、ｈ（ｊ）＜ｇ（ｉ）であれば、ｆ（ｉ）＞ｆ（ｊ）である。これに対して、ｈ（ｉ）＜ｇ（ｊ）であれば、ｆ（ｉ）＜ｆ（ｊ）である。したがって、これら２つの場合には、ｆ（ｉ）＞ｆ（ｊ）であるか否かを判定するために、ｆ（ｉ）とｆ（ｊ）を直接計算しなくてもよい。

処理部１２は、上記のような関係を利用して、ｆ（ｉ）を最小にするｉを決定する。
図１には、処理部１２がプログラム（プログラム）を実行したときの一部の処理（ｆ（ｉ）を最小にするｉで示される状態変数の値を変える部分の処理）の流れの例が示されている。

ステップＳ１：処理部１２は、乱数値であるｒ［１］を生成し、ｉ_ｍｉｎ＝１とする。ｉ_ｍｉｎは、ｆ（ｉ）を最小にするｉの候補を示す識別番号であり、更新候補の状態変数を示す。また、処理部１２は、入力された問題情報（式（１）のＷ_ｉｊ、ｂ_ｉ、ｃ）、ｘ_ｉ（ｉ＝１,２，…,ｎ）の初期値に基づいて、式（２）により、ｘ_１が変化したときのエネルギーの変化量（ΔＥ_１）を計算する。また、処理部１２は、ｇ（ｉ_ｍｉｎ）＝ｇ（１）とｈ（ｉ_ｍｉｎ）＝ｈ（１）を計算する。処理部１２は、ｒ［１］が属する区間のｌｂ［ｋ］、ｕｂ［ｋ］を記憶部１１から取得して（読み出して）、ｇ（１）＝β×ｍａｘ（０，ΔＥ_１）＋ｌｂ［ｋ］、ｈ（１）＝β×ｍａｘ（０，ΔＥ_１）＋ｕｂ［ｋ］を計算する。

ステップＳ２：処理部１２は、ｘ_ｉ（ｉ＝２,３,…,ｎ）を選択する。ｎは状態変数の数である。なお、ステップＳ１の処理の直後は、ｉ＝２である。
ステップＳ３：処理部１２は、ｒ［ｉ］を生成する。

ステップＳ４：処理部１２は、ｒ［ｉ］が属する区間のｌｂ［ｋ］とｕｂ［ｋ］を記憶部１１から取得する（読み出す）。１×Ｄ≦ｒ［ｉ］≦２×Ｄ－１の場合、図１の例のように、ｋ＝２の区間のｌｂ［２］とｕｂ［２］が取得される。

ステップＳ５：処理部１２は、選択したｘ_ｉの値が変化した場合のエネルギーの変化量であるΔＥ_ｉを計算するとともに、ｇ（ｉ）＝β×ｍａｘ（０，ΔＥ_ｉ）＋ｌｂ［ｋ］、ｈ（ｉ）＝β×ｍａｘ（０，ΔＥ_ｉ）＋ｕｂ［ｋ］を計算する。

ステップＳ６：処理部１２は、ｈ（ｉ）とｇ（ｉ_ｍｉｎ）との比較結果、または、ｇ（ｉ）とｈ（ｉ_ｍｉｎ）の比較結果に基づいて、選択したｘ_ｉを更新候補とするか否かを決定する。すなわち、ｉをｉ_ｍｉｎとするか、これまでのｉ_ｍｉｎを維持するかが決定される。

前述の関係を用いて、ｈ（ｉ）＜ｇ（ｉ_ｍｉｎ）の場合、ｆ（ｉ_ｍｉｎ）＞ｆ（ｉ）と判定される。このため、処理部１２は、ｉをｉ_ｍｉｎとする。一方、ｈ（ｉ_ｍｉｎ）＜ｇ（ｉ）の場合、ｆ（ｉ_ｍｉｎ）＜ｆ（ｉ）と判定される。このため、処理部１２は、これまでのｉ_ｍｉｎを維持する。

なお、これら以外の場合（ｇ（ｉ）～ｈ（ｉ）の区間と、ｇ（ｉ_ｍｉｎ）～ｈ（ｉ_ｍｉｎ）の区間が重なっている場合）、処理部１２は、ｆ（ｉ_ｍｉｎ）とｆ（ｉ）を計算する。そして、処理部１２は、ｆ（ｉ_ｍｉｎ）＞ｆ（ｉ）であれば、ｉをｉ_ｍｉｎとし、ｆ（ｉ_ｍｉｎ）≦ｆ（ｉ）であれば、これまでのｉ_ｍｉｎを維持する。

ステップＳ７，Ｓ８：処理部１２は、ｉ＝ｉ＋１とし、ｉ≦ｎであるか否かを判定する。ｉ≦ｎであると判定された場合には、ステップＳ２からの処理が繰り返され、ｉ≦ｎではないと判定された場合、ステップＳ９の処理が行われる。

ステップＳ９：処理部１２は、ｉ_ｍｉｎにより識別される更新候補の状態変数であるｘ_ｉｍｉｎを、更新対象の状態変数として決定し、その値を変化（フリップ）する。
処理部１２は、ステップＳ１～Ｓ９の処理を繰り返すことで、エネルギーが極小になるｎ個の状態変数の値を探索する。

処理部１２は、疑似焼き鈍し法を行う場合、所定の温度パラメータ変更スケジュールにしたがって、βの値を大きくしていく（Ｔの値を小さくしていくことに相当する）。そして、処理部１２は、ステップＳ１～Ｓ９の処理が所定回数、繰り返された場合に得られたｎ個の状態変数の値を、組合せ最適化問題の計算結果として出力する（たとえば、図示しない表示装置に表示する）。なお、処理部１２は、ステップＳ９の処理のたびに、式（１）で表される評価関数の値（エネルギー）を更新し、これまでの最小エネルギーとなった場合のエネルギーとｎ個の状態変数の値とを保持しておいてもよい。その場合、処理部１２は、ステップＳ１～Ｓ９の処理が所定回数、繰り返された後に記憶されている最小エネルギーに対応するｎ個の状態変数の値を、計算結果として出力してもよい。

処理部１２は、レプリカ交換法を行う場合、それぞれ異なる温度パラメータの値が設定された複数のレプリカのそれぞれにおいて、図１に示したステップＳ１～Ｓ９の処理を行う。そして、処理部１２は、ステップＳ１～Ｓ９の処理が一定回数繰り返されるごとに、レプリカ交換を行う。たとえば、処理部１２は、複数のレプリカのうち２つをランダムに選択して、選択された２つのレプリカの間で、レプリカ間のエネルギー差や温度パラメータの値の差に基づいた所定の交換確率で、温度パラメータの値を交換する。処理部１２は、たとえば、各レプリカにおいてステップＳ９の処理のたびに、式（１）で表される評価関数の値（エネルギー）を更新し、これまでの最小エネルギーとなった場合のエネルギーとｎ個の状態変数の値とを保持する。処理部１２は、そして、各レプリカにおいてステップＳ１～Ｓ９の処理が所定回数、繰り返された後に記憶されている最小エネルギーのうち、全レプリカにおいて最小のエネルギーに対応する配置情報を、計算結果として出力する。

以上のような第１の実施の形態のデータ処理装置１０及びデータ処理方法によれば、予め記憶したノイズ関数の区間の下限値及び上限値を用いて、各状態変数についての評価値（ｆ（ｉ））がとり得る下限評価値（ｇ（ｉ））と上限評価値（ｈ（ｉ））が計算される。そして、ｆ（ｉ）がこれまでで最小の状態変数である更新候補についてのｇ（ｉ_ｍｉｎ）とｈ（ｉ_ｍｉｎ）を用い、ｈ（ｉ）とｇ（ｉ_ｍｉｎ）、ｇ（ｉ）とｈ（ｉ_ｍｉｎ）の比較結果から、ｘ_ｉを更新候補とするか判定できる場合がある。すなわち、ｆ（ｉ）の計算をする頻度を少なくできる。ｆ（ｉ）は前述のように、対数計算により得られるノイズ関数を含むため、計算に時間がかかるが、第１の実施の形態のデータ処理装置１０及びデータ処理方法によれば、ｆ（ｉ）の計算頻度を少なくできるため、組合せ最適化問題の計算時間を短縮できる。

（第２の実施の形態）
図３は、第２の実施の形態のデータ処理装置のハードウェア例を示すブロック図である。

データ処理装置２０は、たとえば、コンピュータであり、ＣＰＵ２１、ＲＡＭ２２、ＨＤＤ２３、ＧＰＵ（Graphics Processing Unit）２４、入力インタフェース２５、媒体リーダ２６及び通信インタフェース２７を有する。上記ユニットは、バスに接続されている。

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

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

ＨＤＤ２３は、ＯＳ（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 (7)

1. 複数の状態変数を含むイジング型の評価関数の値が極小または極大となる前記複数の状態変数の値を、乱数値を用いて探索する処理をコンピュータに実行させるプログラムであって、
   前記乱数値の範囲を分割した複数の区間のそれぞれについて、前記乱数値に基づく対数計算により得られるノイズ関数の下限値と上限値を記憶した記憶部から、前記評価関数に含まれる前記複数の状態変数のそれぞれについて生成した前記乱数値が属する区間の前記下限値と前記上限値を、前記複数の状態変数のそれぞれについて取得し、
   前記複数の状態変数のうちの１つの状態変数の値が変化した場合の前記評価関数の値の変化量と０のうちの大きい方の値と温度パラメータとの積と、前記状態変数について取得した前記下限値との和である下限評価値と、前記積と前記状態変数について取得した前記上限値との和である上限評価値とを、前記複数の状態変数のそれぞれについて計算し、
   前記下限評価値と前記上限評価値に基づいて、前記複数の状態変数のそれぞれの間の、前記積と前記ノイズ関数の和による評価値の大小関係を判定することで、前記評価値が最小となる更新対象の状態変数を決定し、
   決定した前記更新対象の状態変数の値を変える、
   処理を繰り返すことで、前記評価関数の値が極小または極大になる前記複数の状態変数の値を探索する、
   処理をコンピュータに実行させるプログラム。
2. 前記複数の状態変数のうちの第１の状態変数について計算した前記下限評価値が、前記複数の状態変数のうちの第２の状態変数について計算した前記上限評価値よりも大きい場合、前記第２の状態変数の前記評価値は、前記第１の状態変数の前記評価値よりも小さいと判定し、
   前記第１の状態変数について計算した前記上限評価値が、前記第２の状態変数について計算した前記下限評価値よりも小さい場合、前記第１の状態変数の前記評価値は、前記第２の状態変数の前記評価値よりも小さいと判定する、
   処理をコンピュータに実行させる請求項１に記載のプログラム。
3. 前記第１の状態変数について計算した前記下限評価値が、前記第２の状態変数について計算した前記上限評価値以下であり、前記第１の状態変数について計算した前記上限評価値が、前記第２の状態変数について計算した前記下限評価値以上である場合、前記第１の状態変数と前記第２の状態変数のそれぞれについて、前記評価値を計算し、
   前記第１の状態変数について計算した前記評価値と、前記第２の状態変数について計算した前記評価値との比較結果に基づいて、前記大小関係を判定する、
   処理を前記コンピュータに実行させる請求項２に記載のプログラム。
4. 前記乱数値の範囲を分割した前記複数の区間のそれぞれについて、前記乱数値に基づく対数計算により得られる前記ノイズ関数の前記下限値と前記上限値とを計算し、前記記憶部に記憶する、
   処理を前記コンピュータに実行させる請求項１乃至３の何れか一項に記載のプログラム。
5. 所定の繰り返し回数ごとに、前記温度パラメータの値を変化させる、処理を前記コンピュータに実行させる請求項１乃至４の何れか一項に記載のプログラム。
6. 複数の状態変数を含むイジング型の評価関数の値が極小または極大となる前記複数の状態変数の値を、乱数値を用いて探索するデータ処理方法であって、
   コンピュータが、
   前記乱数値の範囲を分割した複数の区間のそれぞれについて、前記乱数値に基づく対数計算により得られるノイズ関数の下限値と上限値を記憶した記憶部から、前記評価関数に含まれる前記複数の状態変数のそれぞれについて生成した前記乱数値が属する区間の前記下限値と前記上限値を、前記複数の状態変数のそれぞれについて取得し、
   前記複数の状態変数のうちの１つの状態変数の値が変化した場合の前記評価関数の値の変化量と０のうちの大きい方の値と温度パラメータとの積と、前記状態変数について取得した前記下限値との和である下限評価値と、前記積と前記状態変数について取得した前記上限値との和である上限評価値とを、前記複数の状態変数のそれぞれについて計算し、
   前記下限評価値と前記上限評価値に基づいて、前記複数の状態変数のそれぞれの間の、前記積と前記ノイズ関数の和による評価値の大小関係を判定することで、前記評価値が最小となる更新対象の状態変数を決定し、
   決定した前記更新対象の状態変数の値を変える、
   処理を繰り返すことで、前記評価関数の値が極小または極大になる前記複数の状態変数の値を探索する、
   データ処理方法。
7. 複数の状態変数を含むイジング型の評価関数の値が極小または極大となる前記複数の状態変数の値を、乱数値を用いて探索するデータ処理装置であって、
   前記乱数値に基づく対数計算により得られるノイズ関数の下限値と上限値を、前記乱数値の範囲を分割した複数の区間のそれぞれについて記憶した記憶部と、
   前記記憶部から、前記評価関数に含まれる前記複数の状態変数のそれぞれについて生成した前記乱数値が属する区間の前記下限値と前記上限値を、前記複数の状態変数のそれぞれについて取得し、前記複数の状態変数のうちの１つの状態変数の値が変化した場合の前記評価関数の値の変化量と０のうちの大きい方の値と温度パラメータとの積と、前記状態変数について取得した前記下限値との和である下限評価値と、前記積と前記状態変数について取得した前記上限値との和である上限評価値とを、前記複数の状態変数のそれぞれについて計算し、前記下限評価値と前記上限評価値に基づいて、前記複数の状態変数のそれぞれの間の、前記積と前記ノイズ関数の和による評価値の大小関係を判定することで、前記評価値が最小となる更新対象の状態変数を決定し、決定した前記更新対象の状態変数の値を変える、処理を繰り返すことで、前記評価関数の値が極小または極大になる前記複数の状態変数の値を探索する処理部と、
   を有するデータ処理装置。

Images