JP6851460B2 - 最適解判定方法、最適解判定プログラム、非一時的記録媒体及び最適解判定装置 - Google Patents

Description

本発明は最適解判定方法、最適解判定プログラム及び最適解判定装置に係り、特に組合せ最適化問題における解の最適性を判定する技術に関する。

近年、ビッグデータに対するデータマイニングなど、データ解析のニーズが高まってきている。重要なデータ解析分野のひとつとして、組合せ最適化問題がある（例えば古くから知られている巡回セールスマン問題などが含まれる）。

組合せ最適化問題は、ＮＰ（Non-deterministic Polynomial time）完全あるいはＮＰ困難な難しい問題が多く含まれる。即ち、一般的に問題の規模が大きくなると計算量が指数関数以上のオーダーで爆発するので、網羅的な全探索による解決がほとんど不可能である。

そこで、ヒューリスティックに近似解を求める手法が多く開発され、適用されている。しかし、全探索をしていないので、得られた解が「局所解」に留まり、真の最適解かどうかを判定できないという問題がある。

この問題に直接応えるため、解に何らかのスコアを付与し、統計解析を併用することで、解の最適性を見極めようという研究があり、例えば、特に極値統計解析によって最大スコアを推定する手法などがある（非特許文献１）。

しかしながら、場合によっては複雑な制約条件が課せられる組合せ最適化問題においては、良質な統計解析の前提となる解空間から解候補を大量に一様抽出する手段の提供が難しく、実用化には至っていないと考えられる（非特許文献２）。

また、組合せ最適化問題に関連して、ＺＤＤ（Zero-suppressed Binary Decision Diagram）と呼ばれるデータ構造があり、フロンティア法と称される構築アルゴリズムによって、非常に大規模な組合せ集合を効率的に列挙索引化できることが分かり、近年盛んに研究されている（非特許文献３）。ＺＤＤでは、組合せ最適化問題の解候補全体を描出し、そこから効率的に解候補を一様抽出できる。

また、創薬分野では、近年発達したＮＧＳ（Next Generation Sequencer：次世代シーケンサ）などにより、大量の遺伝子データ、例えばＲＮＡ（ribonucleic acid）発現行列データが取得できるようになった。そうして得られたビッグデータの解析は、バイオインフォマティクスとして注目されている。例えば、生体機能を踏まえた薬剤の作用機序などを解明しようという試みがある。そのひとつに、遺伝子制御ネットワークの推定がある。遺伝子制御ネットワークとは、遺伝子が相互に発現量を調節するシステムを、ベイジアンネットワークなどの確率的グラフモデルとして捉える解析手法である。

特許文献１には、遺伝子発現の時系列研究から遺伝子間のネットワーク関係を推定するためのノンパラメトリック回帰によるベイジアンモデルを使用する技術が開示されている。

国際公開第２００４／０４７０２０号

Golden, B.L. and Alt, F.B. "Interval estimation of a global optimum for large combinatorial problems", Naval Research Logistics Quarterly, 26, 69-77 (1979) Giddings, A.P., Rardin, R.L., and Uzsoy, R. "Statistical optimum estimation techniques for combinatorial optimization problems: a review and critique", Journal of Heuristics, 20(3), 329-358 (2014) Iwashita, H., Kawahara, J., and Minato, S. "ZDD-Based Computation of the Number of Paths in a Graph", TCS Technical Report, TCS-TR-A-12-60, Hokkaido University, September 18, 2012.

遺伝子制御ネットワークでは、遺伝子をノード、制御関係をエッジとするグラフを考える。そして、グラフ構造によって得られているＲＮＡ発現行列データをどの程度説明し得るかを計算して、得られているデータに最も適合するグラフを探索しようという、グラフ探索（グラフマイニング）問題を解きたい。しかしながら、遺伝子数Ｎに対して、可能なグラフ数は2^(N^2)通りあり、Ｎの増大につれて超指数関数的に発散する。また、ベイジアンネットワークモデルでは、非循環有向グラフ制約（DAG制約：Directed acyclic graph制約）という複雑な制約条件も考えなければならない。遺伝子数がある程度大きくなると、ＺＤＤでも解空間全体を描出するには困難を伴う。

しかも、創薬プロジェクトは、人体に供される可能性がある薬剤を生み出そうというもので、結果の判明まで長い時間と莫大なコストを要し、産業上特に、解析結果の妥当性に対する関心が高い。よって、取得したＲＮＡ発現行列データを巧く説明するグラフ構造の最適性が判定できれば、産業上非常に有用である。

特許文献１に記載の遺伝子制御ネットワークを構築する方法は、いくつかの遺伝子に対して遺伝子発現の経時変化データを取得し、ベイジアン推定法の修正と、発現遺伝子間の原因及び結果の関係を判断する。ベイジアン推定法の修正は、経時変化データを用いて発現遺伝子間の因果関係を推定することを含む。特許文献１に記載の方法は、発現の遅い遺伝子の変化は、発現の早い遺伝子の変化の原因となる可能性があまりないという仮定に基づいてベイジアン推定及びノンパラメトリック回帰を修正して信頼性の高いネットワーク解を提供するものであるが、上記仮定が全ての遺伝子制御ネットワークに適用できる保証はない。

本発明はこのような事情に鑑みてなされたもので、組合せ最適化問題における解の最適性の判定を効率的かつ精度よく行うことができる最適解判定方法、最適解判定プログラム及び最適解判定装置を提供することを目的とする。

上記目的を達成するために一の態様に係る発明は、組合せ最適化問題における解の最適性をコンピュータにより判定する最適解判定方法であって、組合せ最適化問題における解の最適性をコンピュータにより判定する最適解判定方法であって、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出する第１のステップと、第１のステップにより一様抽出した第１の複数の解のそれぞれに対応する第１の複数の評価値を取得する第２のステップと、取得した第１の複数の評価値に基づいて、第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第３のステップと、解空間に属する解のうち少なくとも１つの解を解候補として取得する第４のステップと、解候補に対応する評価値を取得する第５のステップと、第５のステップにより取得した解候補に対応する評価値が第１の最大評価値の信頼区間内に入るか否かを判定する第６のステップと、第６のステップにおいて解候補に対応する評価値が第１の最大評価値の信頼区間内に入ると判定された場合に、解候補を第１の最適解とする第７のステップと、を含む。

本発明の一の態様によれば、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出し、一様抽出した第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を第１の最大評価値として推定する。そして、解空間に属する解のうち少なくとも１つの解を解候補として取得すると、その解候補に対応する評価値を取得し、取得した評価値が第１の最大評価値の信頼区間内に入るか否か（解候補が第１の最適解の１つか否か）を判定する。これにより、組合せ最適化問題の解の統計的な最適性判定が可能になり、かつ解が十分条件を満たすか否かの最適性判定が可能である。

本発明の他の態様に係る最適解判定方法において、第１の複数の解は、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解であり、第３のステップは、Ｕ×Ｖ個の解をＶ個のブロックに分け、ブロック毎にＵ個の解の評価値の区分最大値をＶ個取得し、Ｖ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして第１の最大評価値を推定することが好ましい。ブロック毎のＵ個の解の評価値の区分最大値をＶ個取得し、Ｖ個の区分最大値を用いて、区分最大値が一般極値分布（ＧＥＶ：generalized extreme value distribution）に従うものとして最大評価値（第１の最大評価値）を最尤推定する。

本発明の更に他の態様に係る最適解判定方法において、第７のステップによる判定結果を出力する第８のステップを更に含むことが好ましい。

本発明の更に他の態様に係る最適解判定方法において、演算コストは小さいが解の精度が低い第１の探索法と、第１の探索法よりも演算コストは大きいが解の精度が高い第２の探索法とを有し、第４のステップは、最初に第１の探索法により探索された第１の解候補を入力し、第１の解候補の評価値が第１の最大評価値の信頼区間内に入らない場合のみ、第２の探索法により探索された第２の解候補を入力することが好ましい。

最初に演算コストは小さいが解の精度が低い第１の探索法により探索された第１の解候補を入力し、その第１の解候補が十分条件を満たさず十分性判定に失敗した場合、第１の探索法によるヒューリスティック探索が不十分なために、他の最適解が存在することが示唆されたことになる。その場合、演算コストは大きいが解の精度が高い第２の探索法による第２の解候補の探索に切り替え、第２の探索法により探索される第２の解候補の十分性を判定する。

本発明の更に他の態様に係る最適解判定方法において、解空間上の解であって、第１の最適解からの解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出する第９のステップと、第２の複数の解のそれぞれに対応する第２の複数の評価値を取得する第１０のステップと、第２の複数の評価値に基づいて、第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第１１のステップと、第１の最適解の評価値を取得する第１２のステップと、第１の最適解の評価値が、第２の最大評価値を超えているか否かを判定する第１３のステップと、を含むことが好ましい。

本発明の更に他の態様によれば、第１の最適解からの解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出し、一様抽出した第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を第２の最大評価値として推定する。そして、第１の最適解の評価値が第２の最大評価値を超えているか否かを判定する。これにより、組合せ最適化問題の解が必要条件を満たすか否かの最適性判定が可能であり、第１の最適解の評価値が第２の最大評価値を超えている場合には、第１の最適解は必要十分条件を満たし、解空間上で最適な解であり、同等解も他には存在しないことになる。

本発明の更に他の態様に係る最適解判定方法において、第２の複数の解は、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解であり、第１１のステップは、Ｐ×Ｑ個の解をＱ個のブロックに分け、ブロック毎のＰ個の解の評価値の区分最大値をＱ個取得し、Ｑ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして第２の最大評価値を推定することが好ましい。ブロック毎のＰ個の解の評価値の区分最大値をＱ個取得し、Ｑ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして最大評価値（第２の最大評価値）を最尤推定する。

本発明の更に他の態様に係る最適解判定方法において、第１３のステップによる判定結果を出力する第１４のステップを更に含むことが好ましい。

本発明の更に他の態様に係る最適解判定方法において、第１３のステップにより第１の最適解の評価値が第２の最大評価値を超えていないと判定されると、一定範囲を拡大し、拡大した一定範囲外の複数の解を第３の複数の解として一様抽出する第１５のステップと、第３の複数の解のそれぞれに対応する第３の複数の評価値を取得する第１６のステップと、第３の複数の評価値に基づいて、第３の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第３の最大評価値とする第１７のステップと、第１の最適解の評価値が、第３の最大評価値を超えているか否かを判定する第１８のステップと、を更に含むことが好ましい。

最適解の１つとして探索された第１の最適解の評価値が第２の最大評価値を超えていない場合、同等の最適解が存在する可能性があるが、この場合には、第１の最適解から離れた解空間（拡大した一定範囲外の解空間）から一様抽出した第３の複数の解に基づいて第３の最大評価値を再度推定する。これにより、探索された最適解から離れた解空間には、同等の最適解が存在しないことを確認することができる。

本発明の更に他の態様に係る最適解判定方法において、第１３のステップにより第１の最適解の評価値が第２の最大評価値を超えていないと判定されると、解空間上の解であって、かつ第１の最適解から一定距離離れた第４の解候補を取得する第１９のステップと、解空間上において、第１の最適解及び第４の解候補のそれぞれからの解空間上の距離が一定範囲外の第４の複数の解を一様抽出する第２０のステップと、第４の複数の解のそれぞれに対応する第４の複数の評価値を取得する第２１のステップと、第４の複数の評価値に基づいて、第４の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第４の最大評価値とする第２２のステップと、第１の最適解の評価値が、第４の最大評価値を超えているか否かを判定する第２３のステップと、を更に含むことが好ましい。

最適解の１つとして探索された第１の最適解の評価値が第２の最大評価値を超えていない場合、第１の最適解と、第１の最適解から一定距離離れた第４の解候補との両方から、それぞれ解空間上の距離が一定範囲外の第４の複数の解を一様抽出し、一様抽出した第４の複数の解に基づいて第４の最大評価値を再度推定する。これにより、探索された第１の最適解及び一定範囲外の解空間の第４の解候補からそれぞれ離れた解空間には、同等の最適解が存在しないことを確認することができる。

本発明の更に他の態様に係る最適解判定方法において、解空間は、第１の制約条件における第１の解空間と第２の制約条件における第２の解空間とを含み、第１の解空間に属する第１の解候補を取得する第２４のステップと、第２の解空間上の複数の解を第５の複数の解として一様抽出する第２５のステップと、第５の複数の解のそれぞれに対応する第５の複数の評価値を取得する第２６のステップと、取得した第５の複数の評価値に基づいて、第５の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第５の最大評価値とする第２７のステップと、第２の解空間における解候補であって、第２４のステップで取得した第１の解候補からの解空間上の距離が近い解を近傍解として取得する第２８のステップと、近傍解の評価値を取得する第２９のステップと、近傍解の評価値が第５の最大評価値の信頼区間内に入るか否かを判定する第３０のステップと、を含むことが好ましい。

本発明の更に他の態様によれば、第１の制約条件における第１の解空間における解候補の最適性を判定する場合、その解候補からの解空間上の距離が近い近傍解であって、第２の制約条件における第２の解空間における近傍解を取得し、近傍解の最適性を判定する。即ち、近傍解の最適性判定を代用し、近傍解が最適であれば、その近傍解に近い距離の第１の解空間における解候補も最適であると推定する。

本発明の更に他の態様に係る最適解判定方法において、第１のステップは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、組合せ可能パターンを縮約して列挙索引化するデータ構造を用い、解空間上の解の総数を算出し、算出した総数以下の乱数を発生させ、発生させた乱数により特定されるパターンに対応する解を抽出することが好ましい。

これにより、非常に大規模で、効率的な手段を用いてもなお全描出が難しい組合せ最適化問題であっても解を一様抽出することができ、解の最適性判定を行うことができる。

本発明の更に他の態様に係る最適解判定方法において、組合せ最適化問題は、遺伝子制御ネットワークの組合せ最適化問題であることが好ましい。

本発明の更に他の態様に係る最適解判定プログラムは、上記の最適解判定方法をコンピュータに実行させる。

更に他の態様に係る発明は、組合せ最適化問題における解の最適性を判定する最適解判定装置であって、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出する解抽出部と、一様抽出した第１の複数の解のそれぞれに対応する第１の複数の評価値を取得する第１の評価値取得部と、取得した第１の複数の評価値に基づいて、第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第１の最大評価値推定部と、解空間に属する解のうち少なくとも１つの解を解候補として取得する解取得部と、解候補に対応する評価値を第１の評価値取得部から取得し、取得した解候補に対応する評価値が第１の最大評価値の信頼区間内に入るか否かを判定し、解候補に対応する評価値が第１の最大評価値の信頼区間内に入ると判定された場合に、解候補を第１の最適解とする第１の判定部と、を備える。

本発明の更に他の態様に係る最適解判定装置において、解空間上の解であって、第１の最適解からの解空間上の距離が一定範囲外の複数の解を第２の複数の解として一様抽出する第２の解抽出部と、第２の複数の解のそれぞれに対応する第２の複数の評価値を取得する第２の評価値取得部と、第２の複数の評価値に基づいて、第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第２の最大評価値推定部と、第１の最適解に対応する評価値を第２の評価値取得部から取得し、取得した第１の最適解に対応する評価値が、第２の最大評価値を超えているか否かを判定する第２の判定部と、を更に備えることが好ましい。

本発明によれば、組合せ最適化問題の解の統計的な最適性判定が可能になり、組合せ最適化問題における解の最適性の判定を効率的かつ精度よく行うことができる。

ＲＮＡ発現行列データを示す図表及び遺伝子制御ネットワークを示す図 本発明の特徴を示す概念図 本発明に係る最適解判定装置のハードウェア構成を示すブロック図 本発明の第１の実施形態を示す機能ブロック図 要素「Ａ，Ｂ，Ｃ，Ｄ」からなる全体集合「Ｇ_∧」を示す図 ４つの部分集合（Ｇ_（１），Ｇ_（２），Ｇ_（３），Ｇ_（４））を示す図 ３つの部分集合（Ｇ_（１），Ｇ_（３），Ｇ_（４））の選択に対応する「パターン」を示す図 部分集合の選択に対応する「パターン」の３つの例と、各「パターン」が条件を満たしているか否かの判定結果とを示す図 図５から図８に示した集合分割問題における全てのパターン及び判定結果を、「２分グラフ」で網羅的に表現した図 フロンティア法の「枝刈り」により組合せ可能パターンが縮約される様子を示す図 フロンティア法の「節共有」により組合せ可能パターンが縮約される様子を示す図 フロンティア法の「枝刈り」及び「節共有」により組合せ可能パターンを縮約した結果を示す図 ノード＝｛Ａ，Ｂ，Ｃ｝で、エッジの数が３本のグラフ集合をＺＤＤ表現した図及び特定のパス（Ａ⇔Ｂ⇒Ｃ）に対応するグラフを示す図 ＺＤＤ表現したグラフ集合の全体のグラフ数（採用総数）の「数え上げ」を示す図 ＺＤＤ表現したグラフ集合から任意の指定番号のグラフを取り出す方法を示す図 本発明の第２の実施形態を示す機能ブロック図 ある遺伝子制御ネットワーク推定に対して本発明を適用した例を示すグラフ 本発明に係る最適解判定方法の第１の実施形態を示すフローチャート 本発明に係る最適解判定方法の第２の実施形態を示すフローチャートであり、特に最適十分性判定に成功した後に行われる最適必要性判定の処理を示すフローチャート 本発明に係る最適解判定方法の第４の実施形態の要部を示すフローチャート 本発明に係る最適解判定方法の第５の実施形態を示すフローチャート 制約条件Ｃにおける第１の解空間と制約条件Ｃ_＋における第２の解空間とを模式的に表す図

以下、添付図面に従って本発明に係る最適解判定方法、最適解判定プログラム及び最適解判定装置の好ましい実施の形態について説明する。

＜本発明の概要＞
組合せ最適化問題における解を探索する方法として、創薬分野に応用可能な遺伝子制御ネットワークを例に説明する。遺伝子制御ネットワークとは、遺伝子間の協調関係を有向グラフとして表現することで、例えば薬剤の作用機序などを読み解くための応用などが期待されている。

まず、基本的な前提条件について説明する。

(1) 図１は、ＲＮＡ発現行列データを示す図表及び遺伝子制御ネットワークを示す。

図１において、Ａ、Ｂ、…、Ｚは遺伝子であり、Ｘ１、Ｘ２、…、Ｘｎはサンプルであり、遺伝子の数とサンプルの数との積だけデータが存在する。このＲＮＡ発現行列データを取得する。ＲＮＡ発現行列データは、ＮＧＳによるカバレッジデータであってもよいし、マイクロアレイによるシグナルデータでもよい。

ＲＮＡ発現行列では、Ｍ個の細胞株等に対して、Ｎ個の遺伝子のＲＮＡ発現量が計測されており、データx(m,n)は、細胞株ｍにおける遺伝子ｎの発現量を示す。したがって、ＲＮＡ発現行列ＤはＭ×Ｎの数値行列データである。

(2) 複数の遺伝子の関連性は、遺伝子制御ネットワークとして表すことができる。以下、遺伝子制御ネットワークをグラフＧとして表す。図１に示すようにグラフＧは、エッジ（矢印で示すノード（遺伝子）間の制御関係）の集合であり、例えばg1={(A,B),(A,C),(C,D)}は、「遺伝子AからB」「遺伝子AからC」「遺伝子CからD」の３つの制御関係が存在することを示す。グラフＧは、多数のサンプルでのＲＮＡ発現行列データＤから推定可能である。

ただし、グラフＧには問題に応じて何らかの制約Ｃが課せられる。制約はモデル上、もしくは事前知識によって課せられる。例えばベイジアンネットワークモデルでは循環グラフを表現できないため、グラフＧは循環グラフであってはならない。（つまり、例えば「(A,B),(B,A)」や「(A,B),(B,C),(C,A)」といった部分集合を含んではいけない。）また、事前知識によってスケールフリーネットワーク性（ノードの次数分布がべき乗則に適合すること）が期待されていれば、そのような制約を設けることも考えられる。

(3) 評価関数Ｓ（Ｄ，Ｇ）を用意する。これは、グラフＧがデータＤをどれだけ説明できているかを定量化したものである。例えば、前述のg1では、「遺伝子AからB」の制御関係に対しては、「x(m,B)=F(x(m,A))が当てはまるかどうか」を定量化する。Fは制御関係のモデル関数であり、定量化は、例えば罰則付き最大尤度（AIC（Akaike's Information Criterion）やBIC（Bayesian information criterion））が使われる。

(4) 最適化問題を解く何らかのヒューリスティックな手法により、最適と思われるグラフＧ_1を獲得し、それに対する評価値Ｓ_1を取得する。

上記は遺伝子制御ネットワークの推定の例であるが、(a) 何らかのデータＤに基づき、(b) 制約Ｃに従う集合Ｇを考え、(c) 評価関数Ｓの最大化（もしくは最小化）を試みて、(d) 特定の集合Ｇ_1を取得する、というところがポイントであって、かつ、組合せ最適化問題では一般的な構造である。

したがって、本発明は、様々な他の組合せ最適化問題に対しても適用し得る。例えば、細胞の排他的かつ被覆的な遺伝子変異を探索する問題がある。これによって、例えば癌にとって重要な遺伝子変異又は作用機序を推定しようという応用が知られている。

(1) Ｍ個の細胞に対して、Ｎ種類のＳＮＰ（Single Nucleotide Polymorphism）における変異データｘを取得する。データx(m,n)は、細胞mにおけるＳＮＰ変異nの有無を示す。ここで、ＳＮＰとはＤＮＡ（deoxyribonucleic acid）のうち変異が入り易い位置のことを言う。

(2) 遺伝子座のセットを集合Ｇとして表す。集合ＧはＳＮＰの集合であり、例えばg1={1,3,4}は、「ＳＮＰ1,3,4」の３つのＳＮＰに注目することを示す。集合Ｇの要素はＮ種類のＳＮＰのいずれかであり、集合Ｇ全体でＳＮＰの組合せパターンに対応する。ただし、制約Ｃとして、集合ＧによってＭ個の細胞株を排他的に被覆することが要求される。即ち、「ＳＮＰ1に変異を持つ細胞の集合M1」「ＳＮＰ3に変異を持つ細胞の集合M3」「ＳＮＰ4に変異を持つ細胞の集合M4」を考えたとき、M1,M3,M4は互いに重複要素となる細胞を共有してはならず、かつ、M1,M3,M4全体ですべての細胞を網羅していなければならない。

(3) 評価関数Ｓ（Ｇ）を用意する。これは、集合Ｇの何らかの特性を定量化したものである。例えば、ＳＮＰの集合Ｇを、予め用意した遺伝子制御ネットワークＦに載せたときの適合性などを評価することが考えられる。

(4) 何らかのヒューリスティックな手法により、最適と思われるグラフＧ_1を獲得し、それに対する評価値Ｓ_1を取得する。

上記(4)としては、例えば、グリーディ・ヒルクライミング法、焼きなまし法、タブーサーチ、遺伝的アルゴリズムなどが知られており、そのいずれを用いても良い。

また、本発明はいわゆるバイオインフォマティクス以外の分野にも適用し得る。例えば遺伝子制御ネットワーク推定はベイジアンネットワークとして一般化されるので、多数の製品の様々な特性を測定してデータ化し、その特性同士の因果関係を推定する手法としても利用できる。組合せ最適化問題としては、例えばナップザック問題や巡回セールスマン問題などが知られ、様々な分野にて応用されており、本発明はそのいずれにも適用できる。

さて、通常のヒューリスティック探索は、上記(4) でアルゴリズムは終了するので、ヒューリスティック探索により獲得したグラフＧ_1が、真の最適値かどうかを判断できなかった。

例えば遺伝子制御ネットワーク問題の場合、グラフＧ_1がＲＮＡ発現行列データＤに対して真に最も適合しているのかどうかを判断できなかった。そのため、多額のコストを要する介入実験に踏み込むためには、例えばグラフＧ_1をバイオロジストが精査して妥当性を判断するなどの属人的で不確実な工程を要していた。

そこで、本発明は、最適と思われるグラフＧ_1の最適性を見極めることができるようにする。

図２は、本発明の特徴を示す概念図である。

ヒューリスティック探索により探索した最適と思われるグラフＧ_1の最適性を見極めるために、ヒューリスティック探索により探索したグラフの抽出個数を超える個数の解を想定した場合の評価値（解空間全体の解の評価値（「スコア」とも言う））のうちの最大評価値（第１の最大評価値）Ｚを推定する。尚、第１の最大評価値Ｚの具体的な推定方法については後述する。

続いて、最適と思われるグラフＧ_1（ローカル解）が、推定した第１の最大評価値Ｚの信頼区間内に入るか否かを判定する。そして、グラフＧ_1が第１の最大評価値Ｚの信頼区間内に入っていれば、グラフＧ_1は、探索空間（解空間）の全域での第１の最適解（グローバル解の一つ）であると判定できる。上記のような判定（解が十分条件を満たすか否かの最適十分性判定）が、本発明の特徴の一つである。

ただし、上記の最適十分性判定の場合、ローカル解がグローバル解か否かの判定は可能であるが、唯一のグローバル解か否かの判定はできない。したがって、グラフＧ_1（ローカル解）の最適十分性判定に成功しても同等解が他にも存在する可能性があり、探索に「未練」が残る。

本発明は、グラフＧ_1の最適十分性判定に成功した場合、そのグラフＧ_1からの解空間上の距離が一定範囲外の解空間（部分空間）上の解の評価値のうちの最大評価値（第２の最大評価値）Ｗを推定する。尚、第２の最大評価値Ｗの具体的な推定方法については後述する。

続いて、グラフＧ_1が、推定した第２の最大評価値Ｗを超えているか否かの判定を更に行う。そして、グラフＧ_1が第２の最大評価値Ｗを超えていれば、部分空間には、グラフＧ_1と同等解が存在しないことになり、そのグラフＧ_1は唯一のグローバル解であると判定できる。上記のような判定（解が必要条件を満たすか否かの最適必要性判定）が、本発明の他の特徴の一つである。

＜最適解判定装置＞
［装置構成］
図３は本発明に係る最適解判定装置のハードウェア構成を示すブロック図である。

図３に示す最適解判定装置１０は、コンピュータによって構成されており、主として各構成要素の動作を制御する中央処理装置（ＣＰＵ：Central Processing Unit）１２と、装置の制御プログラムが格納されたり、プログラム実行時の作業領域となる主メモリ１４と、液晶ディスプレイ、ＣＲＴディスプレイ等のモニタ装置２８の表示を制御するグラフィックボード１６と、ネットワーク５０と接続される通信インターフェース（通信Ｉ／Ｆ）１８と、本発明に係る最適解判定プログラムを含む各種のアプリケーションソフト、及び後述する最適解の判定結果等を保存するハードディスク装置２０と、ＣＤ−ＲＯＭドライブ２２と、キーボード３０のキー操作を検出して指示入力としてＣＰＵ１２に出力するキーボードコントローラ２４と、位置入力装置としてのマウス３２の状態を検出してモニタ装置２８上のマウスポインタの位置やマウス３２の状態等の信号をＣＰＵ１２に出力するマウスコントローラ２６とから構成されている。

また、ネットワーク５０には、ＲＮＡ発現行列データを保存するデータベース４０が接続されている。ＲＮＡ発現行列データは、図１に示したように複数の細胞株（サンプル：Ｘ１、Ｘ２、…、Ｘｎ）における複数の遺伝子（Ａ、Ｂ、…、Ｚ）のＲＮＡ発現量を示す数値行列データである。また、ＲＮＡ発現量は、図示しないＮＧＳ（Next Generation Sequencer）などによりサンプルから取得されたものである。

最適解判定装置１０は、通信インターフェース１８を介してデータベース４０にアクセスし、必要なＲＮＡ発現行列データを取得することができる。尚、ＲＮＡ発現行列データは、外部のデータベース４０に格納されたものを使用する場合に限らず、ＲＮＡ発現行列データをハードディスク装置２０に保存し、ハードディスク装置２０に保存されたＲＮＡ発現行列データを使用するようにしてもよい。

［第１の実施形態］
図４は、図３に示した最適解判定装置１０のＣＰＵ１２の機能を示す機能ブロック図であり、本発明の第１の実施形態を示す機能ブロック図である。

ＣＰＵ１２は、ハードディスク装置２０に格納された最適解判定プログラムを実行することより各種の処理部として機能し、図４に示す第１の実施形態では、解抽出部１００、第１の評価値取得部及び第２の評価値取得部として機能する評価値取得部１０２、解入力部１０４、第１の最大評価値推定部１０６、第１の比較部１０８及び第１の判定部１１０としての機能を有する。

解抽出部１００は、組合せ最適化問題（遺伝子制御ネットワーク）の解空間上の解（グラフ）を一様抽出する部分であり、本例では、ＺＤＤ（Zero-suppressed binary Decision Diagram）を用いたパス列挙索引化アルゴリズムによりグラフＧを列挙索引化する。制約Ｃの下でＺＤＤを構築することで、グラフＧの総数と、集合｛Ｇ｝の任意の要素を一様抽出できるようになる。遺伝子制御ネットワークでは、遺伝子をノード、制御関係をエッジとするグラフを考える。

尚、ここでいう「解」は、「許容解（実行可能解）」を意味し、最適解とは限らないが、実行不可能ではない解のことを指す。すなわち（実行が不可能な）不適な解はあらかじめ排除されている。

また、解抽出部１００は、「枝刈り」及び「節共有」の少なくとも一方を用いて、組合せ最適化問題における組合せ可能なパターンを縮約して列挙索引化するデータ構造とする。ここで、「枝刈り」とは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減する処理をいう。また、フロンティア法の「節共有」とは、組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減する処理をいう。

尚、これはＺＤＤだけに限定されるものではなく、組合せ最適化問題に応じて例えばＢＤＤ（Binary Decision Diagram）、あるいはπＤＤ（Permutation Decision Diagram）などのＺＤＤに類似する変形データ構造を用いても構わない。

また、「枝刈り」の「不適かどうか」の判定には制約Ｃを用いる。例えば遺伝子制御ネットワークの場合では、採用済のエッジによって既に循環が発生すれば、残りのエッジを考慮しなくても不適であることが確定する。また、「節共有」の「共通かどうか」の判定においても、制約Ｃのもとでパターンのうち考慮済の部分と残りの要素とを勘案して、共通化できるかどうかを判断する。例えばエッジの数のみを考慮している場合、採用済のエッジの本数が同じであれば、「節共有」を適用できる。「枝刈り」及び「節共有」のアルゴリズムは、ＺＤＤの場合、フロンティア法として知られているので、それに従えばよい。

尚、「枝刈り」のみを備えた手法として、分岐限定法によって解を列挙する手法を用いてもよい。ＺＤＤのうち「節共有」を無効化することで分岐限定法による列挙も可能となるし、「枝刈り」を無効化することも考えられる。しかし、望ましくはＺＤＤのように「枝刈り」及び「節共有」の両方とも利用することで、それによって効率的に解を列挙できる。

解を一様抽出する手段としては、例えばランダム生成も考えられる。即ち、解候補をランダムに生成（例：グラフのエッジのありなしを乱数で決定）し、Ｇの制約を満たしていなければ再生成を繰り返す。制約条件が簡単であれば、ランダム生成の段階でその制約を織り込んでも構わない。例えばエッジ数を所定数以下にする場合、エッジありを選ぶ数に上限を設けておけばよい。一方、例えば循環グラフを禁止したい場合、単純なランダム生成において制約をかけるのは難しいので、ランダム生成されたグラフが循環しているかどうかを判定することが考えられる。しかし、理論的には、ランダム生成でも十分な実施数を重ねることで統計的推定は可能だが、特に判定及び再生成を繰り返す手法で、ランダム生成がカバーする解空間に対して制約下の解空間が小さい場合、例えば解空間のサイズを１：Ｎとすれば、１個の解を生成するのに平均してＮ個のランダム生成が必要になるので、効率が悪い。

したがって、本手法では十分なサンプルサイズの確保が重要なため、特にＺＤＤの導入は効果が大きいと期待される。

＜ＺＤＤの概説＞
次に、ＺＤＤ及びフロンティア法について具体的に説明する。

まず、組合せ最適化問題の一種である集合分割問題へのＺＤＤの適用を考える。

集合分割問題は、ある全体集合に対する部分集合の列が与えられたとき、その幾つかを選んで、「選んだ部分集合同士に重複がない（相互排他）」、かつ「元の全体集合を尽くす（全体被覆）」のようなパターン（組合せ）を作れるか、という問題である。

集合は、「要素の集まり」として定義される。図５に示すように「Ｇ_∧」は全体集合を示し、「Ａ，Ｂ，Ｃ，Ｄ」が「要素」に対応する。

要素を固定し、要素の有り又は無しを「１又は０」に割り振ると、例えば、図６に示すように部分集合が決まる。

「部分集合を含むかどうか」を符号化すれば、図７に示すように「パターン」を表現できる。図７の例では、３つの部分集合（Ｇ_（１），Ｇ_（３），Ｇ_（４））の選択に対応する「パターン」が示されている。

続いて、パターンごとに「条件（被覆性及び排他性）を満たしているか」を判定する。

図８には、部分集合の選択に対応する「パターン」の３つの例と、各「パターン」が条件を満たしているか否かの判定結果とが示されている。

図８に示すように部分集合（Ｇ_（１），Ｇ_（３），Ｇ_（４））の選択に対応する「パターン」は条件を満たし、部分集合（Ｇ_（２），Ｇ_（３），Ｇ_（４））及び部分集合（Ｇ_（３），Ｇ_（４））の選択に対応する「パターン」は条件を満たさない。

部分集合（Ｇ_（２），Ｇ_（３），Ｇ_（４））は要素「Ｃ」が重複し、排他性を満たさず、部分集合（Ｇ_（３），Ｇ_（４））は要素「Ａ」が不足して被覆性を満たさないからである。

さて、本例の集合分割問題における全てのパターン及び判定結果は、図９に示すように「２分グラフ」で網羅的に表現することができる。

本例の部分集合は４個だったため、全てのパターンは、１６個（＝２^４）であるが、部分集合がＮ個の場合、全てのパターンは、２^Ｎとなる。「２分グラフ」の表現では、「２のべき乗」で枝及び葉が増えるという問題がある。

組合せ問題は、「組合せ爆発」が発生し、有限時間で最適解を探索することができなくなることが知られているが、フロンティア法の「枝刈り」及び「節共有」の縮約技術により、所定の条件を満たす組合せ集合の効率的な（実用的な）「数え上げ」を実現することができる。

図１０は、フロンティア法の「枝刈り」により組合せ可能パターンが縮約される様子を示す図である。

図１０に示すように部分集合「Ｇ_（１）」と「Ｇ_（２）」とが同時に選択されると、要素「Ａ」が重複し、その時点で（後続選択に関わらず）不適が確定する。同様に部分集合「Ｇ_（１）」と「Ｇ_（２）」のどちらも選択しないと、要素「Ａ」を漏らすことになり、不適が確定する。そして、不適が確定すると、その時点で展開を打ち切り、「判定結果＝０」に即繋げる。

このように「枝刈り」は、パターン選択の途中で不適が確定すると、その時点で展開を打ち切ることで組合せ可能パターンの縮約を図る。

図１１は、フロンティア法の「節共有」により組合せ可能パターンが縮約される様子を示す図である。

図１１に示すように２つの部分集合「Ｇ_（１）」と「Ｇ_（３）」とが選択される場合と、１つの部分集合「Ｇ_（２）」が選択される場合とでは、いずれも「要素「Ａ」及び「Ｂ」を１回だけ含む」こと、かつ「要素「Ａ」及び「Ｂ」は、もはや含まないようにする」こと、が条件になるので、後続選択による採否判定が完全に一致する。

そして、この場合には、別々に展開せずに、まとめて扱ってよい、つまり、同じ「節」を共有して構わない。

このように「節共有」は、複数のパターン選択の後続が同一ならば、これらをまとめて扱うことで組合せ可能パターンの縮約を図る。

図１２は、フロンティア法の「枝刈り」及び「節共有」により組合せ可能パターンを縮約した結果を示す図である。

図１２に示すように、１６個のパターン（図９）に拡がるはずの枝及び葉を大幅に削減でき、網羅的に１つずつ判定する場合と完全に一致する判定結果を取得することができる。

図１３は、ノード＝｛Ａ，Ｂ，Ｃ｝で、エッジの数が３本のグラフ集合をＺＤＤ表現した図である。

図１３において、実線で示した「１枝」を通る場合のみ、そのエッジを含むグラフ、（点線で示した「０枝」もしくは飛ばしたエッジは含まない）であって、最終的に「１端」に到達したら、そのグラフを採用する。（「０端」に到達したグラフは採用しない。）
ここで、例えば、（Ａ⇔Ｂ⇒Ｃ）のパスに対応するグラフは、図１３の右側の太線の矢印で示した経路で表される。

図１４は、全体のグラフ数（採用総数）の「数え上げ」を示す図である。

図１４に示すように採用総数は、判定結果を示す「１端」に「１」を付与し、「１端」から最上位のＺＤＤノードまで逆順に辿って数え上げることで算出することができる。

数え上げは、下層ＺＤＤノードから各々の枝先の付与数を加算し、加算した数を自身に付与する。これを最上位のＺＤＤノードまで繰り返し、最上位のＺＤＤノードに付与された数値が、全体の採用総数（「１端」に到達するパスの総数）となる。本例の場合、全体の採用総数は、２０個になる。

このように採用総数を算出する「数え上げ」は、ＺＤＤの重要な性質の一つである。

図１５は、任意の指定番号のグラフを取り出す方法を示す図である。

「数え上げ」後、任意の番号（本例の場合、１〜２０の範囲内の番号）が指定されると、指定された番号（指定番号）にしたがって根から下ることで、指定番号に対応するグラフを取り出すことができる。

例えば、「１２番」のグラフを取り出す場合、図１５の太線の矢印にしたがって最上位のノードから下る。まず、最上位のノード（Ａ，Ｂ）から「０枝」又は「１枝」のうち指定番号を含む枝に進む。本例では、「１枝」側の枝を進み、最上位のノード（Ａ，Ｂ）から下層のノード（Ａ，Ｃ）（図１４上で右側のノード（Ａ，Ｃ））に下る。「１枝」側の枝を進む場合、指定番号から「０枝」側の個数を引く。本例では、「１２番」のグラフを取り出すため、指定番号「１２」から「０枝」側の個数「１０」が引かれ、「２」になる。これを「１端」に到達するまで繰り返すことで、図１５の太線の矢印で示すパス（グラフ）を取り出すことができる。尚、図１５に示す「１２番」のグラフは、（Ａ⇒Ｂ⇔Ｃ）のパスに対応するグラフである。

このように採用総数の「数え上げ」後、採用総数内の任意の番号を指定すると、指定した番号により一意に特定されるグラフを取り出すことができる。これにより採用総数以下の乱数を発生させることで、解空間上の解（グラフ）を「一様抽出」することができる。この解空間上の解の「一様抽出」は、ＺＤＤの重要な性質の一つである。

図４に戻って、解抽出部１００は、ＺＤＤを用いたパス列挙索引化アルゴリズムにより遺伝子制御ネットワークの解空間上の解（グラフＧ）を一様抽出する。解空間上のグラフＧの総数は、ＺＤＤの重要な性質の一つである「数え上げ」により行うことができ、また、「数え上げ」により取得したグラフＧの総数以下の乱数（例えば、Ｍ系列（maximal length sequence）を用いた疑似乱数）を発生させ、各乱数により指定された指定番号に対応するグラフＧを取り出す（グラフＧを一様抽出する）。

評価値取得部（第１の評価値取得部）１０２は、解抽出部１００により抽出されたグラフＧに評価値を付与する。例えば、グラフＧがＲＮＡ発現行列データＤをどれだけ説明できているかを定量化した評価関数Ｓ（Ｄ，Ｇ）を用意しておき、評価値取得部１０２は、抽出されたグラフＧに対応する評価値Ｓを評価関数Ｓ（Ｄ，Ｇ）に基づいて取得し、取得した評価値ＳをグラフＧに付与する。評価関数Ｓ（Ｄ，Ｇ）は、データベース４０に保存されたＲＮＡ発現行列データに基づいて評価値取得部１０２が作成してもよいし、予めＲＮＡ発現行列データに基づいて作成され、例えばデータベース４０に保存された評価関数Ｓ（Ｄ，Ｇ）を使用するようにしてもよい。

解取得部１０４は、遺伝子制御ネットワークの最適化問題を解く何らかのヒューリスティックな手法により獲得された最適と思われる解（解候補）（以下、「グラフＧ_1」という）を取得する。また、入力したグラフＧ_1に対する評価値Ｓ_1は、評価値取得部１０２により取得することができる。尚、ヒューリスティックな手法としては、例えば、グリーディ法、ヒルクライミング法、焼きなまし法、タブーサーチ、遺伝的アルゴリズムなどが知られており、そのいずれを用いても良い。また、ヒューリスティックな探索は、本装置により行ってもよいし、外部の装置により行ってもよく、解取得部１０４は、何らかのヒューリスティックな手法により獲得された解候補（グラフＧ_1）を取得する。

第１の最大評価値推定部１０６は、解空間上の解（グラフＧ）の最大評価値を推定する部分であり、本例では、解抽出部１００により一様抽出された複数の解（第１の複数の解）の評価値に基づいて、第１の複数の解の抽出個数を超える個数の解を想定した場合の評価値のうちの最大評価値（第１の最大評価値Ｚ）を推定する。

具体的には、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解（第１の複数のグラフＧ）を一様抽出し、各々のグラフＧに評価値Ｓを与える。ここで、Ｕはブロックサイズであり、Ｖはブロック個数である。Ｕ、Ｖは、ある大きな数として設定する。例えばＵ、Ｖともに10,000に設定しても良く、この場合、一様抽出されるグラフＧの個数は、１億個（＝10,000×10,000)となる。

第１の最大評価値推定部１０６は、Ｕ×Ｖ個のグラフＧをＶ個のブロックに分け、ブロック毎のＵ個のグラフＧの評価値のうちの区分最大値を取得する。したがって、区分最大値は、Ｖ個取得すことができる。そして、Ｖ個の区分最大値が一般極値分布（ＧＥＶ：generalized extreme value distribution）に従うものとして最大評価値（第１の最大評価値Ｚ）を最尤推定する。

第１の最大評価値は、統計学的な裏付けを伴うものである。解空間上のグラフ｛Ｇ｝は本来は有限集合であり、厳密には縮退してしまうが、グラフＧの総数が十分に大きいため連続分布近似が適用できる。その場合、評価値Ｓに明らかに上限が存在するので、適切なＵ、Ｖの設定によってガンベル型になることが期待され、真の第１の最大評価値Ｚを信頼区間付で推定できる。

第１の比較部１０８は、解取得部１０４が取得した解候補（グラフＧ_1）に対応する評価値Ｓ_1を評価値取得部１０２から取得し、取得した評価値Ｓ_1と推定した信頼区間付の第１の最大評価値Ｚとを比較する。

第１の判定部１１０は、第１の比較部１０８による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、第１の最大評価値Ｚの信頼区間内に入るか否かを判定する。第１の最大評価値Ｚの信頼区間内に収まっていれば、グラフＧ_1は、解空間の全域での第１の最適解の一つ（「グラフＧ_1は十分」）であることが分かる。

仮に、Ｚ≫Ｓ_1であれば、両者の評価値の差を解空間上の距離に変換し、現在推定しているグラフＧ_1が真の最適解（第１の最大評価値Ｚに対応する解）からどのくらい離れているかを推定してもよい。

即ち、ヒューリスティック探索で探索された解候補（グラフＧ_1）が、真の最適値かどうかを判断できるようになる。

［第２の実施形態］
図１６は、図３に示した最適解判定装置１０のＣＰＵ１２の機能を示す機能ブロック図であり、本発明の第２の実施形態を示す機能ブロック図である。尚、図１６において、図４に示した第１の実施形態と共通する部分には同一の符号を付し、その詳細な説明は省略する。

図１６に示す第２の実施形態はて第１の実施形態と比較して第２の最大評価値推定部１１２、第２の比較部１１４及び第２の判定部１１６が主として追加されている。

第１の実施形態では、解候補（Ｇ_1）の最適性十分性判定に成功しても、解候補（Ｇ_1）以外にも同等に確からしい別の解候補が存在する可能性は排除できない。

第２の実施形態は、解候補（Ｇ_1）以外には同等に確からしい別の解候補が存在しないことを判定可能にするものである。

図１６において、第１の実施形態により解候補（グラフＧ_1）の最適性十分性判定に成功すると、その後、解抽出部１００は、グラフＧ_1からの解空間上の距離が一定範囲外の解を列挙索引化するＺＤＤを構築する。このＺＤＤの構築は、ＺＤＤを構築する際のフロンティア法において、構築途上のグラフＧとグラフＧ_1との共通又は非共通エッジをカウントすることで実現できる。解抽出部１００は、再構築したＺＤＤを用いたパス列挙索引化アルゴリズムにより遺伝子制御ネットワークの解空間上の第２の複数の解（グラフＧ）であって、グラフＧ_1からの解空間上の距離が一定範囲外の第２の複数のグラフＧを一様抽出する。

第２の最大評価値推定部１１２は、一様抽出された第２の複数のグラフにそれぞれ対応する評価値に基づいて、第２の複数の解の抽出個数を超える個数の解を想定した場合の評価値のうちの最大評価値（第２の最大評価値Ｗ）を推定する。

具体的には、第１の最大評価値推定部１０６による第１の最大評価値Ｚの推定と同様の手法により第２の最大評価値Ｗを推定する。即ち、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解（第２の複数のグラフＧ）を一様抽出し、各々のグラフＧに評価値Ｓを与える。ここで、Ｐはブロックサイズであり、Ｑはブロック個数である。Ｐ、Ｑは、第１の最大評価値Ｚを推定する際に一様抽出したＵ、Ｖと同じでもよいし、異なっていてもよい。

第２の最大評価値推定部１１２は、Ｐ×Ｑ個のグラフＧをＱ個のブロックに分け、ブロック毎のＵ個のグラフＧの評価値のうちの区分最大値を取得する。したがって、区分最大値は、Ｑ個取得すことができる。そして、Ｑ個の区分最大値が一般極値分布に従うものとして第２の最大評価値Ｗを最尤推定する。

第２の比較部１１４は、解取得部１０４が取得した解候補（グラフＧ_1）に対応する評価値Ｓ_1を評価値取得部（第２の評価値取得部）１０２から取得し、取得した評価値Ｓ_1と推定した信頼区間付の第２の最大評価値Ｗとを比較する。

第２の判定部１１６は、第２の比較部１１４による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、第２の最大評価値Ｗ（信頼区間付の第２の最大評価値Ｗの範囲）を超えているか否かを判定する。

第１の最大評価値Ｚは、解空間全体の最大値を推定したものであり、第２の最大評価値Ｗは、部分空間の最大値を推定したものなので、基本的にはＷ≦Ｚは明らかである（サンプルサイズ等によっては確率的にＷ＞Ｚとなる場合もある）。

その上で、Ｓ_1≫Ｗ（信頼区間を外れる等）であれば、グラフＧ_1は第１の最大評価値を与えるグラフであって、しかもグラフＧ_1から解空間上で離れた範囲には、同等以上にＲＮＡ発現行列データＤを説明できるグラフ構造はないと判断できる。

尚、グラフＧ_1から離間させる一定範囲の距離は、事前に設定する必要があるが、これはグラフの特性から設定してもよいし、経験的に設定してもよいし、例えばＳ_1≫Ｗになるところまで徐々に離れる距離を大きくしていってもよい。例えば、十分大きな距離の設定値から２分検索などによって効率的に適切な距離を繰り返し探索しても構わない。

言うまでもないことだが、距離ゼロであれば、Ｗ≒Ｚとなり、第１の実施形態と変わらない結果になろうし、ゼロ以外の最短の距離の場合はグラフＧ_1のみしか排除しないため、ある程度の距離を設定しないと、結果には大きな違いは出にくいと想定される。

これにより、第２の最大評価値Ｗを超えていれば、グラフＧ_1は、解空間の全域での唯一の最適解（「グラフＧ_1は必要十分」）であることが分かる。即ち、ヒューリスティック探索で得られたグラフＧ_1以外に、最適値に相当するグラフが存在しないと判断できるようになる。

尚、第２の実施形態は、第１の実施形態において、グラフＧ_1の最適十分性の判定に成功した場合に、続けて行うのが通例であるが、何らかの理由によって決められたグラフＧが存在する場合、そのグラフＧから離れた範囲の解候補とグラフＧとを直接比較する手段として利用されてもよい。

従来のヒューリスティック探索では、例えば初期値をランダムに変えるとか、データにノイズを与えるなどの工夫で繰り返し探索を行うなどの方法もあったが、それ自体もヒューリスティックな判定法であるのに対して、本発明は統計的根拠に裏付けられた手法である。

図１７は、ある遺伝子制御ネットワーク推定に対して本発明を適用した例を示すグラフである。

図１７に示すグラフの横軸は乖離度（距離）であり、縦軸は評価値である。点線は、獲得グラフに対する到達評価値（第１の最大評価値Ｚ）である。また、白抜きの線は、各乖離度に対する最適値の推定値である。

尚、ここでいう「乖離度」とは一般に距離またはそれに準じる指標である。例えば、集合をバイナリ配列として表現する場合はハミング距離としてもよい。問題に応じて、編集距離などの他の距離指標、あるいは、準距離、半距離等の指標を用いてもよい。また、それらに近しい乖離度指標を定義してもよい。

乖離度ゼロでの最適値推定範囲は、到達評価値を包含しているため、到達評価値は最適値であると判定された。

乖離度を増やすに連れて推定範囲に対応する値は徐々に低下し、乖離度５〜６で推定範囲（第２の最大評価値Ｗ）が到達評価値を外れたので、乖離度５〜６以上の範囲には獲得グラフ同等以上の評価値を有するグラフは存在しないことと判定された。

また、実線はヒューリスティック探索で得られた解候補に対応する実測値を示すが、本発明が実態を正しく推定できていることを示す。

＜最適解判定方法＞
［第１の実施形態］
図１８は、本発明に係る最適解判定方法の第１の実施形態を示すフローチャートである。

図１８において、図４に示した解抽出部１００は、組合せ最適化問題の一つある遺伝子制御ネットワークの解空間上の解（グラフＧ）を一様抽出する（ステップＳ１０（第１のステップ））。本例では、前述したようにＺＤＤを用いたパス列挙索引化アルゴリズムによりグラフＧを一様抽出する。

続いて、評価値取得部１０２は、一様抽出されたグラフＧに対して評価値Ｓを付与する（ステップＳ１２）。例えば、グラフＧがＲＮＡ発現行列データＤをどれだけ説明できているかを定量化した評価関数Ｓ（Ｄ，Ｇ）を用意しておき、評価値取得部１０２は、一様抽出されたグラフＧに対応する評価値Ｓを評価関数Ｓ（Ｄ，Ｇ）に基づいて取得し、取得した評価値ＳをグラフＧに付与する。

第１の最大評価値推定部１０６は、一様抽出された複数のグラフＧ（第１の複数の解）に対して、ステップＳ１２により付与された第１の複数の評価値Ｓを取得し（ステップＳ１４、第２のステップ）、取得した第１の複数の評価値Ｓに基づいて、第１の複数のグラフＧの抽出個数を超える個数の解を想定した場合の評価値のうちの最大評価値（第１の最大評価値）Ｚを推定する（ステップＳ１６、第３のステップ）。

具体的には、第１の複数のグラフＧは、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個のグラフであり、第１の最大評価値推定部１０６は、Ｕ×Ｖ個のグラフＧをＶ個のブロックに分け、ブロック毎のＵ個のグラフＧの評価値の区分最大値をＶ個取得し、Ｖ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして最大評価値（信頼区分付の第１の最大評価値Ｚ）を最尤推定する。

続いて、解取得部１０４は、遺伝子制御ネットワークの最適化問題を解く何らかのヒューリスティックな手法により獲得された最適と思われる解候補（グラフＧ_1）を取得する（ステップＳ１８、第４のステップ）。ヒューリスティックな手法としては、例えば、グリーディ法、ヒルクライミング法、焼きなまし法、タブーサーチ、遺伝的アルゴリズムなどが知られており、そのいずれを用いても良い。

評価値取得部１０２（第１の評価値取得部）により、取得したグラフＧ_1に対する評価値Ｓ_1を取得する（ステップＳ２０、第５のステップ）。

第１の比較部１０８は、ステップＳ２０で取得した評価値Ｓ_1と、ステップＳ１６で推定した信頼区間付の第１の最大評価値Ｚとを比較する（ステップＳ２２、第６のステップ）。

第１の判定部１１０は、第１の比較部１０８による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、信頼区間付の第１の最大評価値Ｚの範囲に入るか否かを判定する（ステップＳ２４、第７のステップ）。即ち、第１の判定部１１０は、評価値Ｓ_1が信頼区間付の第１の最大評価値の信頼区間内に入っている場合には、グラフＧ_1は最適解としての十分条件を満たしている（最適十分性あり）と判定する。

第１の判定部１１０による判定結果は、図３に示したモニタ装置２８に表示され、又はハードディスク装置２０に保存され、又は図示しないプリンタにプリント出力される（ステップＳ２６、第８のステップ）。尚、第１の判定部１１０による判定結果は、最適十分性の有無に限らず、最適十分性がない場合には、評価値Ｓ_1と第１の最大評価値Ｚとの差を解空間上の距離に変換し、現在推定しているグラフＧ_1が真の最適解（第１の最大評価値Ｚに対応する解）からどのくらい離れているかを判定してもよい。

第１の実施形態によれば、ヒューリスティック探索で探索された解候補（グラフＧ_1）が、真の最適値としての最適十分性を有するか否かを判定することができる。また、最適十分性の判定に失敗した場合でも、ヒューリスティック探索が不十分なために、他の最適解が存在することが示唆されたことになり、更に失敗程度によって、ある程度は最適解に近いことを主張したり、あるいは、最適解は別にあることを留意した上で、解候補を利用したりしても良い。即ち、最適十分性判定の成否に関わらず、最適十分性判定の情報は有用である。

［第1の実施形態の変形例］
本手順の最適十分性判定に失敗した場合、ヒューリスティックな同じ探索を異なる設定などで繰り返しても良いが、予めヒューリスティック探索を行う複数の探索法（第１の探索法、第２の探索法等）を用意し、複数の探索法を切り替えて使用する方法が考えられる。

例えば、ヒューリスティックな探索法として、演算コストは小さい（探索時間は短い）が解の精度が低い第１の探索法と、第１の探索法よりも演算コストは大きい（探索時間は長い）が解の精度が高い第２の探索法とを準備しておき、最初に第１の探索法により探索した第１の解候補（グラフＧ_1）の最適十分性を判定し、最適十分性の判定に失敗した場合のみ、第２の探索法により探索した第２の解候補（グラフＧ_1）の最適十分性を判定する。

これらの探索法の切り替えは、ヒューリスティックな探索法同士で切り替えてもよいし、近似性の保証がある程度ある近似アルゴリズムや厳密解を求める手法などへ切り替えてもよい。また、探索法は３つ以上用意して順次切り替えても構わない。探索法の切り替えは、方法自体に限らず、同一の探索法の収束判定等によって実現しても構わない。例えば、繰り返しサーチによって精度を高める探索法において、所定回の探索結果を第１の最大評価値Ｚで判定し、第１の最大評価値Ｚの信頼区間内に達するまで探索を繰り返しても良い。

また、この場合、先行して最適十分性を判定するための第１の最大評価値Ｚを取得しておいても構わない。

［第２の実施形態］
図１９は、本発明に係る最適解判定方法の第２の実施形態を示すフローチャートであり、特に図１８に示した第１の実施形態による最適十分性判定に成功した後に行われる最適必要性判定の処理に関して示している。

図１９において、図１６に示した解抽出部１００（第２の解抽出部）は、遺伝子制御ネットワークの解空間上の第２の複数のグラフＧであって、最適十分性判定に成功したグラフＧ_1からの解空間上の距離が一定範囲外の複数（第２の複数）のグラフＧを一様抽出する（ステップＳ３０、第９のステップ）。尚、第２の複数のグラフＧの一様抽出は、グラフＧ_1からの解空間上の距離が一定範囲外の解を列挙索引化するＺＤＤを構築し、構築したＺＤＤを用いたパス列挙索引化アルゴリズムにより行うことができる。

続いて、評価値取得部１０２は、一様抽出された第２の複数のグラフＧに対してそれぞれ評価値Ｓを取得する（ステップＳ３２、第１０のステップ）。グラフＧに対する評価値Ｓの取得は、図１８に示した第１の実施形態のステップＳ１２と同様に行うことができる。

次に、第２の最大評価値推定部１１２は、ステップＳ３２により取得した第２の複数の評価値Ｓに基づいて、第２の複数のグラフＧの抽出個数を超える個数の解を想定した場合の評価値のうちの最大評価値（第２の最大評価値）Ｗを推定する（ステップＳ３４、第１１のステップ）。

具体的には、第２の複数のグラフＧは、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個のグラフであり、第２の最大評価値推定部１１２は、Ｐ×Ｑ個のグラフＧをＱ個のブロックに分け、ブロック毎のＰ個のグラフＧの評価値の区分最大値をＱ個取得し、Ｑ個の区分最大値を用いて、区分最大値が一般極値分布に従うものとして最大評価値（信頼区分付の第２の最大評価値Ｗ）を最尤推定する。

続いて、第２の比較部１１４は、最適十分性判定に成功したグラフＧ_1の評価値Ｓ_1と、ステップＳ３４で推定した信頼区間付の第２の最大評価値Ｗとを比較する（ステップＳ３６、第１２のステップ）。

第２の判定部１１６は、第２の比較部１１４による比較結果に基づいてグラフＧ_1の評価値Ｓ_1が、信頼区間付の第２の最大評価値Ｗの範囲を超えているか否かを判定する（ステップＳ３８、第１３のステップ）。即ち、第２の判定部１１６は、評価値Ｓ_1が信頼区間付の第２の最大評価値の範囲を超えている場合には、グラフＧ_1と同等解は他には存在せず、グラフＧ_1は唯一の最適解としての必要条件を満たしている（最適必要性あり）と判定する。

第２の判定部１１６による判定結果は、図３に示したモニタ装置２８に表示され、又はハードディスク装置２０に保存され、又は図示しないプリンタにプリント出力される（ステップＳ４０、第１４のステップ）。

第２の実施形態によれば、ヒューリスティック探索で探索された解候補（グラフＧ_1）が、真の最適値としての最適十分性判定に成功すると、更にグラフＧ_1の最適必要性を判定するため、探索されたグラフＧ_1から離れた解空間には、同等の評価値を有する解が存在しないことを確認することができる。

［第３の実施形態］
本発明に係る最適解判定方法の第３の実施形態は、図１９に示した第２の実施形態において、最適十分性判定に失敗した場合の処理を含む。

即ち、最適十分性判定に失敗した場合、図１９に示したステップＳ３０（第１５のステップ）において、グラフＧ_1からの一定範囲を解空間上で拡大し、拡大した一定範囲外の解空間上で第３の複数のグラフＧ_3（第３の複数の解）を一様抽出する。

そして、拡大した一定範囲外の解空間上で抽出した第３の複数のグラフＧに対応する第３の複数の評価値の取得（ステップＳ３２、第１６のステップ）、第３の複数の評価値に基づいて第３の複数のグラフＧの個数を超える個数の解を想定した場合の第３の最大評価値Ｗの推定（ステップＳ３４、第１７のステップ）、グラフＧ_1の評価値Ｓ_1と第３の最大評価値Ｗとの比較（ステップＳ３６）、及びグラフＧ_1の最適必要性の判定（ステップＳ３８、第１８のステップ）等を再度実行する。

最適十分性判定に失敗した場合には、一定範囲を拡大させる上記の処理を最適十分性判定に成功するまで、一定範囲を徐々に拡大して複数回繰り返してもよい。

［第４の実施形態］
本発明に係る最適解判定方法の第４の実施形態は、図１９に示した第２の実施形態において、最適十分性判定に失敗した場合の他の処理を含む。

最適十分性判定に失敗した場合（即ち、グラフＧ_1の評価値Ｓ_1が第２の最大評価値Ｗを超えていないとの判定結果が出力された場合）、図２０に示すようにグラフＧ_1からの解空間上の距離が一定距離離れた第４の解候補（グラフＧ_2)であって、第２の最大評価値Ｗに対応する評価値Ｓ_2を有するグラフＧ_2を取得する(ステップＳ４２、第１９のステップ)。

グラフＧ_2の取得は、グラフＧ_1からの距離が一定範囲外の解空間の範囲で解候補（グラフＧ_2)のヒューリスティック探索を実行し、グラフＧ_2の評価値Ｓ_2が、Ｗ≒Ｚと同等となるグラフＧ_2を取得することにより行う。

次に、図１９に示したステップＳ３０の代わりに、グラフＧ_1及びグラフＧ_2のそれぞれからの解空間上の距離が一定範囲外の第４の複数の解を一様抽出する（ステップＳ４４、第２０のステップ）。

その後、図１９に示したステップＳ３２に遷移させ、上記のようにして拡大した一定範囲外の解空間上で抽出した第４の複数のグラフＧに対応する第４の複数の評価値の取得（ステップＳ３２、第２１のステップ）、第４の複数の評価値に基づいて第４の複数のグラフＧの個数を超える個数の解を想定した場合の第４の最大評価値Ｗの推定（ステップＳ３４、第２２のステップ）、グラフＧ_1の評価値Ｓ_1と第４の最大評価値Ｗとの比較（ステップＳ３６）、及びグラフＧ_1の最適必要性の判定（ステップＳ３８、第２３のステップ）等を再度実行する。

そして、最適十分性判定に再度失敗した場合には、新たにグラフＧ_3等を追加し、新たな第４の最大評価値Ｗを推定して最適必要性の判定を繰り返し行う。

尚、一定範囲（意味のある距離）は、組合せ最適化問題のもとの要請から決まる。例えば遺伝子制御ネットワークの場合、距離ｄは異なるエッジの本数を意味するので、作用機序解明等で許容されるエッジの間違い数によって解釈すればよい。エッジの総数がＮ本程度と見込まれる場合、ｄ／Ｎで誤答率を表す指標も考えられるので、例えば誤答率５％を許容してＮ=100が見込まれる場合は、ｄ=5等と設定すればよい。

［第５の実施形態］
解空間全体から解を一様抽出することができない場合がある。解空間全体のＺＤＤが構築できない事例等、そもそも一様抽出手段が確保できない場合である。

その場合に適用可能な本発明に係る最適解判定方法の第５の実施形態について説明する。

図２１は本発明に係る最適解判定方法の第５の実施形態を示すフローチャートである。

まず、制約条件ＣだけではＺＤＤを構築できない場合、ＺＤＤを構築できる、より厳しい制約条件Ｃ_＋を考える。例えばグラフ問題では、制約条件Ｃが「非循環」のみであれば、「全域林」を付与すること等が考えられる。組合せ最適化問題の解空間は、上記の制約条件Ｃ（第１の制約条件）における第１の解空間と、制約条件Ｃ_＋（第２の制約条件）における第２の解空間とを含むものとする。

図２１において、第１の解空間に属する解(制約条件Ｃにおけるグラフ｛Ｇ｝）の中から第１の解候補（グラフＧ_1）を取得する（ステップＳ５０、第２４のステップ）。第１の解空間では、ＺＤＤ等を構築できていなくても、ヒューリスティックな手法等でグラフＧ_1の探索は可能である。

続いて、第２の解空間上の解（制約条件Ｃ_＋におけるグラフ｛Ｇ_＋｝）を列挙し、ステップＳ５０で取得したグラフＧ_1の、制約条件Ｃ_＋における近傍解Ｇ_1_＋を探索する（ステップＳ５２、第２８のステップ）。例えば、制約条件Ｃ_＋におけるＺＤＤは記述できることから、制約条件Ｃ_＋に対してさらに「Ｇ_1からの解空間上の距離が一定範囲内であること」を付け加えてＺＤＤを構築し、そのうち最小距離のものを選ぶことなどにより近傍解Ｇ_1_＋を探索できる。

次に、近傍解Ｇ_1_＋の評価値Ｓ_1_＋を取得する（ステップＳ５４、第２９のステップ）。近傍解Ｇ_1_＋に対する評価値Ｓ_1_＋の取得は、図１８に示した第１の実施形態のステップＳ１２と同様に行うことができる。

ステップＳ５４で取得した近傍解Ｇ_1 ＋の評価値Ｓ_1 ＋と、信頼区間付の第５の最大評価値Ｚとを比較する（ステップＳ５６）。尚、第５の最大評価値Ｚは、第２の解空間の制約条件Ｃ ＋におけるグラフ｛Ｇ ＋｝の一様抽出に基づいて前述した手法によって推定することできる。即ち、第２の解空間上の複数の解を第５の複数の解（制約条件Ｃ ＋におけるグラフ｛Ｇ ＋｝）として一様抽出し（第２５のステップ）、第５の複数の解のそれぞれに対応する第５の複数の評価値を取得し（第２６のステップ）、取得した第５の複数の評価値に基づいて、第５の複数の解の個数を超える個数の解を想定した場合の第５の最大評価値を推定する（第２７のステップ）。

そして、近傍解Ｇ_1_＋の評価値Ｓ_1_＋と、信頼区間付の第１の最大評価値Ｚとの比較結果に基づいて、グラフＧ_1の最適十分性を判定する（ステップＳ５８、第３０のステップ）。即ち、評価値Ｓ_1_＋が第１の最大評価値Ｚの信頼区間内に入るか否かを判定し、評価値Ｓ_1_＋が第１の最大評価値Ｚの信頼区間内に入り、近傍解Ｇ_1_＋の最適十分性があると判定されれば、グラフＧ_1も最適十分性があると判定する。近傍解Ｇ_1_＋が最適なら、少なくとも「より制約が厳しい条件を課した場合の最適解から、最も近いところに解を見つけている」ことは担保できるからである。

図２２は、制約条件Ｃにおける第１の解空間と、制約条件Ｃ_＋における第２の解空間とを模式的に表す図である。ＺＤＤを構築できない制約条件Ｃにおける第１の解空間上でヒューリスティックな手法等で探索したグラフＧ_1の最適性判定を行う場合、ＺＤＤを構築できる、より厳しい制約条件Ｃ_＋における第２の解空間上で近傍解Ｇ_1_＋（グラフＧ_1からの解空間上の距離が一定範囲内の最小距離のもの）を探索し、この近傍解Ｇ_1_＋の最適性判定を代用し、グラフＧ_1の最適性判定を行う。

第２の解空間は、制約が追加されるので「Ｃ_＋⊆Ｃ」となる（図２２）。なぜなら、制約条件Ｃに含まれる制約は、モデルの前提条件であったり、事前知識による妥当な制約であったりするため、基本的には外すべきではないからである。言うまでもなく、追加される制約は、事前知識等によって問題に対してある程度の妥当性が想定されることが望ましい（が、必ずしも完全な妥当性が担保されていなくても止むを得ない）。

しかし、もし制約を外した方が探索容易で、かつ、制約を外しても妥当性の致命的な毀損には至らない場合であれば、「Ｃ⊆Ｃ_＋」になるような制約、即ち制約の解除を考えても構わない。

いずれにせよ、制約条件ＣとＣ_＋との差及びグラフＧ_1と近傍解Ｇ_1_＋との最小距離ｄは、代用した近傍解Ｇ_1_＋と実際に利用する到達解との乖離度を示す。近傍解Ｇ_1_＋の探索については、例えばヒューリスティックな手法によって、到達解の最小距離ｄの近傍については探索していると言えるのであれば、実用上は問題にならないと考えられる。一方、制約条件ＣとＣ_＋との違いの影響は、それによって評価値がどの程度大きく変化するかに依存する。評価値にそれほど大きな影響を与えないような制約条件で、最小距離dが妥当な範囲であれば、近傍解Ｇ_1_＋の最適性判定を代用し、グラフＧ_1の最適性判定を行う信頼度は高いと考えてよい。以上から、より広い組合せ最適化問題に対して本発明の効果を及ぼすことができるようになる。

［その他］
本実施形態の最適解判定装置１０は、例示に過ぎず、他の構成に対しても本発明を適用することが可能である。各機能構成は、任意のハードウェア、ソフトウェア、或いは両者の組合せによって適宜実現可能である。例えば、上述の最適解判定装置１０の各部における処理をコンピュータに実行させる最適解判定プログラム、そのような最適解判定プログラムを記録したコンピュータ読み取り可能な記録媒体（非一時的記録媒体）に対しても、本発明を適用することが可能である。

また、本実施形態において、例えば、解抽出部１００、評価値取得部１０２、解取得部１０４、第１の最大評価値推定部１０６、第１の比較部１０８、及び第１の判定部１１０等の各種の処理を実行する処理部（processing unit）のハードウェア的な構造は、次に示すような各種のプロセッサ（processor）である。各種のプロセッサには、ソフトウェア（プログラム）を実行して各種の処理部として機能する汎用的なプロセッサであるＣＰＵ（Central Processing Unit）、ＦＰＧＡ（Field Programmable Gate Array）などの製造後に回路構成を変更可能なプロセッサであるプログラマブルロジックデバイス（Programmable Logic Device：ＰＬＤ）、ＡＳＩＣ（Application Specific Integrated Circuit）などの特定の処理を実行させるために専用に設計された回路構成を有するプロセッサである専用電気回路などが含まれる。

１つの処理部は、これら各種のプロセッサのうちの１つで構成されていてもよいし、同種または異種の２つ以上のプロセッサ（例えば、複数のＦＰＧＡ、あるいはＣＰＵとＦＰＧＡの組み合わせ）で構成されてもよい。また、複数の処理部を１つのプロセッサで構成してもよい。複数の処理部を１つのプロセッサで構成する例としては、第１に、クライアントやサーバなどのコンピュータに代表されるように、１つ以上のＣＰＵとソフトウェアの組合せで１つのプロセッサを構成し、このプロセッサが複数の処理部として機能する形態がある。第２に、システムオンチップ（System On Chip：ＳｏＣ）などに代表されるように、複数の処理部を含むシステム全体の機能を１つのＩＣ（Integrated Circuit）チップで実現するプロセッサを使用する形態がある。このように、各種の処理部は、ハードウェア的な構造として、上記各種のプロセッサを１つ以上用いて構成される。

さらに、これらの各種のプロセッサのハードウェア的な構造は、より具体的には、半導体素子などの回路素子を組み合わせた電気回路（circuitry）である。

また、本発明は、プロセッサを有する最適解判定装置であって、プロセッサが、組合せ最適化問題の解空間上の複数の解を第１の複数の解として一様抽出し、一様抽出した第１の複数の解のそれぞれに対応する第１の複数の評価値を取得し、取得した第１の複数の評価値に基づいて、第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を第１の最大評価値として推定し、解空間に属する解のうち少なくとも１つの解を解候補として取得し、取得した解候補に対応する評価値が第１の最大評価値の信頼区間内に入るか否かを判定する最適解判定装置を含む。

更に、本発明は上述した実施形態に限定されず、本発明の精神を逸脱しない範囲で種々の変形が可能であることは言うまでもない。

１０ 最適解判定装置
１２ ＣＰＵ
１４ 主メモリ
１６ グラフィックボード
１８ 通信インターフェース
２０ ハードディスク装置
２２ ＣＤ−ＲＯＭドライブ
２４ キーボードコントローラ
２６ マウスコントローラ
２８ モニタ装置
３０ キーボード
３２ マウス
４０ データベース
５０ ネットワーク
１００ 解抽出部
１０２ 評価値取得部
１０４ 解入力部
１０６ 第１の最大評価値推定部
１０８ 第１の比較部
１１０ 第１の判定部
１１２ 第２の最大評価値推定部
１１４ 第２の比較部
１１６ 第２の判定部
Ｚ 第１の最大評価値
Ｗ 第２の最大評価値

Claims (14)

1. 組合せ最適化問題における解の最適性をコンピュータが以下の各ステップの処理を行うことにより判定する最適解判定方法であって、
   前記組合せ最適化問題の解空間上の複数の解を、組合せ可能パターンを縮約して列挙索引化するデータ構造を用いて、又はランダム生成により第１の複数の解として一様抽出する第１のステップと、
   前記第１のステップにより一様抽出した前記第１の複数の解のそれぞれに対応する第１の複数の評価値を、評価関数に基づいて取得する第２のステップと、
   取得した前記第１の複数の評価値に基づいて、前記第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第３のステップと、
   前記解空間に属する解のうち少なくとも１つの解を、ヒューリスティックな探索法により解候補として取得する第４のステップと、
   前記解候補に対応する評価値を、前記評価関数に基づいて取得する第５のステップと、
   前記第５のステップにより取得した前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入るか否かを判定する第６のステップと、
   前記第６のステップにおいて前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入ると判定された場合に、前記解候補を第１の最適解とする第７のステップと、を含み、
   前記第１の複数の解は、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解であり、
   前記第３のステップは、前記Ｕ×Ｖ個の解をＶ個のブロックに分け、前記ブロック毎にＵ個の解の評価値の区分最大値をＶ個取得し、前記Ｖ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第１の最大評価値を推定する、
   最適解判定方法。
2. 前記第７のステップによる判定結果を出力する第８のステップを更に含む、
   請求項１に記載の最適解判定方法。
3. 演算コストは小さいが解の精度が低い第１の探索法と、前記第１の探索法よりも演算コストは大きいが解の精度が高い第２の探索法とを有し、
   前記第４のステップは、最初に前記第１の探索法により探索された第１の解候補を入力し、前記第１の解候補の評価値が前記第１の最大評価値の範囲に入らない場合のみ、前記第２の探索法により探索された第２の解候補を入力する、
   請求項１又は２に記載の最適解判定方法。
4. 前記解空間上の解であって、前記第１の最適解からの前記解空間上の距離が一定範囲外の複数の解を、組合せ可能パターンを縮約して列挙索引化するデータ構造を用いて、又はランダム生成により第２の複数の解として一様抽出する第９のステップと、
   前記第２の複数の解のそれぞれに対応する第２の複数の評価値を、前記評価関数に基づいて取得する第１０のステップと、
   前記第２の複数の評価値に基づいて、前記第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第１１のステップと、
   前記第１の最適解の評価値を取得する第１２のステップと、
   前記第１の最適解の評価値が、前記第２の最大評価値を超えているか否かを判定する第１３のステップと、を含み、
   前記第２の複数の解は、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解であり、
   前記第１１のステップは、前記Ｐ×Ｑ個の解をＱ個のブロックに分け、前記ブロック毎にＰ個の解の評価値の区分最大値をＱ個取得し、前記Ｑ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第２の最大評価値を推定する、
   請求項１から３のいずれか１項に記載の最適解判定方法。
5. 前記第１３のステップによる判定結果を出力する第１４のステップを更に含む、
   請求項４に記載の最適解判定方法。
6. 前記第１３のステップにより前記第１の最適解の評価値が前記第２の最大評価値を超えていないと判定されると、前記一定範囲を拡大し、前記拡大した一定範囲外の複数の解を第３の複数の解として一様抽出する第１５のステップと、
   前記第３の複数の解のそれぞれに対応する第３の複数の評価値を取得する第１６のステップと、
   前記第３の複数の評価値に基づいて、前記第３の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第３の最大評価値とする第１７のステップと、
   前記第１の最適解の評価値が、前記第３の最大評価値を超えているか否かを判定する第１８のステップと、
   を更に含む、
   請求項４又は５に記載の最適解判定方法。
7. 前記第１３のステップにより前記第１の最適解の評価値が前記第２の最大評価値を超えていないと判定されると、前記解空間上の解であって、かつ前記第１の最適解から一定距離離れた第４の解候補を取得する第１９のステップと、
   前記解空間上において、前記第１の最適解及び前記第４の解候補のそれぞれからの前記解空間上の距離が一定範囲外の第４の複数の解を一様抽出する第２０のステップと、
   前記第４の複数の解のそれぞれに対応する第４の複数の評価値を取得する第２１のステップと、
   前記第４の複数の評価値に基づいて、前記第４の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第４の最大評価値とする第２２のステップと、
   前記第１の最適解の評価値が、前記第４の最大評価値を超えているか否かを判定する第２３のステップと、
   を更に含む、
   請求項４又は５に記載の最適解判定方法。
8. 前記解空間は、第１の制約条件における第１の解空間と第２の制約条件における第２の解空間とを含み、
   前記第１の解空間に属する第３の解候補を取得する第２４のステップと、
   前記第２の解空間上の複数の解を第５の複数の解として一様抽出する第２５のステップと、
   前記第５の複数の解のそれぞれに対応する第５の複数の評価値を取得する第２６のステ
   ップと、
   取得した前記第５の複数の評価値に基づいて、前記第５の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第５の最大評価値とする第２７のステップと、
   前記第２の解空間における解候補であって、前記第２４のステップで取得した前記第３の解候補からの前記解空間上の距離が近い解を近傍解として取得する第２８のステップと、
   前記近傍解の評価値を取得する第２９のステップと、
   前記近傍解の評価値が前記第５の最大評価値の信頼区間内に入るか否かを判定する第３０のステップと、
   を含む請求項１から７のいずれか１項に記載の最適解判定方法。
9. 前記第１のステップは、前記組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び前記組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、前記組合せ可能パターンを縮約して列挙索引化するデータ構造を用い、前記解空間上の解の総数を算出し、前記算出した総数以下の乱数を発生させ、発生させた乱数により特定されるパターンに対応する解を抽出する、
   請求項１から８のいずれか１項に記載の最適解判定方法。
10. 組合せ最適化問題は、遺伝子制御ネットワークの組合せ最適化問題である請求項１から９のいずれか１項に記載の最適解判定方法。
11. 請求項１から１０のいずれか１項に記載の最適解判定方法をコンピュータに実行させる最適解判定プログラム。
12. 請求項１１に記載の最適解判定プログラムを記録したコンピュータ読み取り可能な非一時的記録媒体。
13. 組合せ最適化問題における解の最適性を判定する最適解判定装置であって、
    前記組合せ最適化問題の解空間上の複数の解を、組合せ可能パターンを縮約して列挙索引化するデータ構造を用いて、又はランダム生成により第１の複数の解として一様抽出する解抽出部と、
    前記一様抽出した前記第１の複数の解のそれぞれに対応する評価値を、評価関数に基づいて取得する第１の評価値取得部と、
    前記取得した前記第１の複数の評価値に基づいて、前記第１の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第１の最大評価値とする第１の最大評価値推定部と、
    前記解空間に属する解のうち少なくとも１つの解を、ヒューリスティックな探索法により解候補として取得する解取得部と、
    前記解候補に対応する評価値を、前記評価関数に基づいて前記第１の評価値取得部から取得し、前記取得した前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入るか否かを判定し、前記解候補に対応する評価値が前記第１の最大評価値の信頼区間内に入ると判定された場合に、前記解候補を第１の最適解とする第１の判定部と、を備え、
    前記第１の複数の解は、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解であり、
    前記第１の最大評価値推定部は、前記Ｕ×Ｖ個の解をＶ個のブロックに分け、前記ブロック毎にＵ個の解の評価値の区分最大値をＶ個取得し、前記Ｖ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第１の最大評価値を推定する、
    最適解判定装置。
14. 前記解空間上の解であって、前記第１の最適解からの前記解空間上の距離が一定範囲外の複数の解を、組合せ可能パターンを縮約して列挙索引化するデータ構造を用いて、又はランダム生成により第２の複数の解として一様抽出する第２の解抽出部と、
    前記第２の複数の解のそれぞれに対応する第２の複数の評価値を、前記評価関数に基づいて取得する第２の評価値取得部と、
    前記第２の複数の評価値に基づいて、前記第２の複数の解の個数を超える個数の解を想定した場合の最大評価値を推定し、第２の最大評価値とする第２の最大評価値推定部と、
    前記第１の最適解に対応する評価値を前記第２の評価値取得部から取得し、前記取得した前記第１の最適解に対応する評価値が、前記第２の最大評価値を超えているか否かを判定する第２の判定部と、を更に備え、
    前記第２の複数の解は、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解であり、
    前記第２の最大評価値推定部は、前記Ｐ×Ｑ個の解をＱ個のブロックに分け、前記ブロック毎にＰ個の解の評価値の区分最大値をＱ個取得し、前記Ｑ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第２の最大評価値を推定する、
    請求項１３に記載の最適解判定装置。

Images