JP6751376B2 - 最適解探索方法、最適解探索プログラム及び最適解探索装置 - Google Patents

Description

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

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

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

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

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

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

組合せ問題に対して、ヒューリスティックに近似解を求める手法が多く開発され、適用されている。その中でもよく用いられる手法に、グリーディ・ヒルクライミング法がある。これは、解空間{Ｓ}に対して何らかの初期解S_0を与え、次いで初期解S_0から乖離度１の解候補群{S_0+}を総当たり的に評価して適合度{P_0+｝を定量化し、その中から最大（または最小）な解S_1を選ぶ。これを、適合度が改善しなくなるまで繰り返すという手法である。一方、この手法は局所解に陥りやすい欠点があると言われていて、様々な修正法が提案されている。

そのような改良法のひとつに、次のステップに進む際に、乖離度１ではなく、乖離度Ｎの解候補を列挙する手法がある。しかしながら、この手法では乖離度Ｎが大きくなるに連れて、次の解候補群の数が爆発的に増えるため、特に、複雑な制約条件では、列挙自体が難しいという課題があった。また、解候補群の個数が非常に大きい場合は、それを均等に探索するのが難しいという課題があった。

例えば、特許文献１では、最適となる組合せ状態を求めるために、初期の組合せ状態から出発し、隣接状態と定義された組合せ状態の中から遷移すべき組合せ状態を決定して、順次、組合せ状態の遷移を繰り返してネットワーク構成の探索を行う最適解を求めるデータ処理方法であって、問題固有の状態間の距離である問題固有距離（乖離度）を定義し、探索による評価関数値の改善率が一定値以下となった時点で乖離度の大きな遷移を数回行い、その後、所定回数の探索において乖離度の小さな状態変化に限定して探索を行い、これを繰り返すことによって探索を継続するデータ処理方法が提案されている。

また、組み合わせ最適化問題の探索アルゴリズムにＺＤＤを活用する研究についても、発展途上である。

例えば、特許文献２では、ＺＤＤまたはＢＤＤ（Binary Decision Diagram）（二分決定グラフ）を用いてナップザック（詰め込み）問題を効率的に解く手法が提案されている。

特許文献２に記載の詰め込み支援装置は、詰め込み可能な詰め込み物の数の組合せの制約から生成された、詰め込み可能な詰め込み物の数の組合せの制約を展開した、詰め込み物の数の組合せパターンの集合を表す二分決定グラフを二分決定グラフ情報記憶部に記憶させ、ユーザにより指定された詰め込み物の数の組合せパターンを取得すると、二分決定グラフ情報記憶部に記憶された二分決定グラフを用いた探索により、取得した詰め込み物の数の組合せパターンが詰め込み可能であるか否かを判定する。

特開２０１０−１８６４２５号公報 特許第５９８７５３０号公報

組合せ最適化問題における最適解の探索において、次の解候補を、元の解候補から乖離度の大きな解候補群に対して効率的に探索することが難しく、特に乖離度の大きな解候補群が評価可能な個数を超えた場合には、均等な解候補の探索ができないという問題がある。

特許文献１には、組合せ状態と乖離度（問題固有距離）、それに改善率とを考慮し、遷移すべき組合せ状態を決定することが提案されているが、ある乖離度の隣接状態にある解をどのように把握するかについての効率的な解決手段については検討が尽くされていない。

また、特許文献２では、ＺＤＤまたはＢＤＤを用いて詰め込み問題を効率的に解く手法を提案しているが、直接的に解を描出するための利用であり、ナップザック問題以外の解全体を描出し切れない場合には対応できないという問題がある。

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

上記目的を達成するために一の態様に係る発明は、組合せ最適化問題における最適解をコンピュータにより探索する最適解探索方法であって、組合せ最適化問題の解空間に属する解のうち少なくとも１つの解を第１の解候補として取得する第１のステップと、第１の解候補に評価値を付与する第２のステップと、第１の解候補からの乖離度が第１の範囲以下に収まる解候補群を二分決定グラフとして列挙索引化する第３のステップであって、二分決定グラフは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、組合せ可能パターンを縮約して列挙索引化するデータ構造を有する、第３のステップと、列挙索引化された解候補群から解候補群の一部又は全部を第２の解候補として一様抽出する第４のステップと、抽出した第２の解候補に評価値を付与する第５のステップと、第１の解候補の評価値及び第２の解候補の評価値の１つ以上の評価値に基づいて第１の最適解の探索の終了の是非を判断する第６のステップと、を含み、第１の最適解の探索が終了していないと判断された場合は、第２の解候補の中から選択された１つ以上の解候補であって、第１の解候補と異なる解候補を第１の解候補として更新し、第３のステップから第６のステップの処理を繰り返し、第１の最適解の探索が終了したと判断された場合は、終了と判断された評価値が付与された第１の解候補を、第１の最適解として出力する。

本発明の一の態様によれば、組合せ最適化問題の解空間に属する解のうち少なくとも１つの解を第１の解候補とし、この第１の解候補に評価値を付与するとともに、第１の解候補からの乖離度が第１の範囲以下に収まる解候補群を二分決定グラフとして列挙索引化する。二分決定グラフは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減する（いわゆる「枝刈り」により組合せ可能パターンを縮約する）ステップ、及び組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減する（いわゆる「節共有」により組合せ可能パターンを縮約する）ステップのうちの少なくとも一方を用いて、組合せ可能パターンを縮約して列挙索引化するデータ構造を有し、例えば、ＺＤＤ（Zero-suppressed Binary Decision Diagram）、またはＺＤＤに類似するデータ構造を有する。

そして、列挙索引化された解候補群から解候補群の一部又は全部を第２の解候補として一様抽出し、一様抽出した第２の解候補に評価値を付与し、第１の解候補の評価値及び第２の解候補の評価値の１つ以上の評価値に基づいて第１の最適解の探索の終了の是非を判断する。第１の最適解の探索が終了していないと判断された場合は、第２の解候補の中から選択された１つ以上の解候補であって、第１の解候補と異なる解候補を第１の解候補として更新することで、新たな第１の解候補からの乖離度が第１の範囲以下に収まる解候補群を二分決定グラフとして列挙索引化し、上記の第１の最適解の探索を繰り返す。一方、第１の最適解の探索が終了したと判断された場合は、終了と判断された評価値が付与された第１の解候補を、第１の最適解（最適解）として出力する。

これにより、広い解空間（第１の解候補から乖離度の大きな解候補群）であっても、効率的に解候補群を列挙索引化して最適解を検索することができ、解候補群が評価可能な個数を超えたとしても一様抽出することで均等な探索が可能である。

本発明の他の態様に係る最適解探索方法において、１つ以上の解候補の制約条件を受け付ける第７のステップを含み、第３のステップは、第１の解候補からの乖離度が第１の範囲以下に収まり、かつ制約条件を満たす解候補群を二分決定グラフとして列挙索引化することが好ましい。制約条件を課した中で解探索する場合、ＺＤＤ又はＺＤＤに類似するデータ構造の導入は更に大きな効果が期待できる。

本発明の更に他の態様に係る最適解探索方法において、第１の範囲の乖離度は、１以上二分決定グラフとして列挙索引化が可能な最大の乖離度以下であることが好ましい。

本発明の更に他の態様に係る最適解探索方法において、第１の範囲の乖離度は、一定値又は第１の解候補が更新される毎に変化する値であることが好ましい。

本発明の更に他の態様に係る最適解探索方法において、第６のステップは、第１の解候補の評価値が全ての第２の解候補の評価値以上の場合、第２の解候補と第１の解候補との差が規定値以下の場合、又は第３のステップから第６のステップの処理の繰り返し回数が一定回数に達する場合を、第１の最適解の探索の終了と判断することが好ましい。現在の第１の解候補よりも評価値の高い解候補（第２の解候補）が探索されなくなり、あるいは評価値の改善度合いが小さくなり、あるいは探索し尽くしたからです。

本発明の更に他の態様に係る最適解探索方法において、第１の解候補として更新される第１の解候補と異なる解候補は、第２の解候補の中で最大の評価値が付与された第２の解候補であることが好ましい。これにより、第１の解候補は、現在の第１の解候補よりも評価値の高い解候補に次々と更新され、最終的に第１の最適解に到達することができる。

本発明の更に他の態様に係る最適解探索方法において、第６のステップは、第５のステップで付与した第２の解候補の評価値に基づいて、第２の解候補の個数を超える個数の解を想定した場合の最大の評価値を、第１の最大評価値として推定する第８のステップと、第２の解候補の評価値が第１の最大評価値の信頼区間内に入るか否かを判定する第９のステップと、を含み、第１の解候補の評価値及び第１の最大評価値の信頼区間内に入ると判定された第２の解候補の評価値に基づいて第１の最適解の探索の終了の是非を判断することが好ましい。

本発明の更に他の態様によれば、第５のステップで付与した第２の解候補の評価値に基づいて、第２の解候補の個数を超える個数の解を想定した場合の最大の評価値を、第１の最大評価値として推定し、一様抽出される第２の解候補の評価値が第１の最大評価値の信頼区間内に入るか否か（第２の解候補が解候補群のうちの最大の評価値を有する解候補の１つか否か）を判定する。これにより、組合せ最適化問題の解の統計的な最適性判定が可能になり、かつ一様抽出される第２の解候補が、第１の解候補からの乖離度が第１の範囲以下に収まる解候補群の中で、最大の評価値を有する解候補としての十分条件を満たすか否かの最適性判定が可能である。

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

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

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

本発明の更に他の態様に係る最適解探索方法において、第１の解候補からの乖離度が第１の範囲以下に収まる解候補群に含まれる第２の解候補群であって、第１の最大評価値の信頼区間内に評価値が入ると判定された第２の解候補からの乖離度が、第１の範囲よりも狭い第２の範囲外の第２の解候補群を二分決定グラフとして列挙索引化する第１０のステップであって、二分決定グラフは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、組合せ可能パターンを縮約して列挙索引化するデータ構造を有する、第１０のステップと、列挙索引化された第２の解候補群から第２の解候補群の一部又は全部を第３の解候補として一様抽出する第１１のステップと、抽出した第３の解候補に評価値を付与する第１２のステップと、第１２のステップで付与した第３の解候補の評価値に基づいて、第３の解候補の個数を超える個数の解を想定した場合の最大の評価値を、第２の最大評価値として推定する第１３のステップと、第１の最大評価値の信頼区間内に入った第２の解候補の評価値が、第２の最大評価値を超えているか否かを判定する第１４のステップと、を含むことが好ましい。

本発明の更に他の態様によれば、第１の解候補からの乖離度が第１の範囲以下に収まる解候補群に含まれる第２の解候補群であって、第１の最大評価値の信頼区間内に評価値が入ると判定された第２の解候補からの乖離度が、第１の範囲よりも狭い第２の範囲外の第２の解候補群を二分決定グラフとして列挙索引化する。この二分決定グラフは、ＺＤＤ又はこれに類似するデータ構造を有する。そして、列挙索引化された第２の解候補群から第２の解候補群の一部又は全部を第３の解候補として一様抽出し、第３の解候補に評価値を付与し、第３の解候補の評価値に基づいて、第３の解候補の個数を超える個数の解を想定した場合の最大の評価値を、第２の最大評価値として推定する。第１の最大評価値の信頼区間内に入った第２の解候補の評価値が、第２の最大評価値を超えているか否かを判定することで、組合せ最適化問題の解が必要条件を満たすか否かの最適性判定を可能にしている。第１の最大評価値の信頼区間内に入った第２の解候補の評価値が第２の最大評価値を超えている場合には、その第２の解候補は必要十分条件を満たし、第１の解候補からの乖離度が第１の範囲内では最適な解であり、同等解も他には存在しないことになる。

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

本発明の更に他の態様に係る最適解探索方法において、第１４のステップにより第１の最大評価値の信頼区間内に入った第２の解候補の評価値が、第２の最大評価値を超えていないと判定されると、第２の範囲に代えて第２の範囲を拡大した第３の範囲を適用して第１０のステップから第１４のステップの処理を行うことが好ましい。

第１の解候補からの乖離度が第１の範囲内で最適な解の１つとして探索された第２の解候補の評価値が第２の最大評価値を超えていない場合、第１の解候補からの乖離度が第１の範囲内に同等の解が存在する可能性があるが、この場合には、第２の範囲に代えて第２の範囲を拡大した第３の範囲を適用し、第１の解候補からより離れた解候補群を使用して第２の最大評価値を再度推定する。これにより、最適な解の１つとして探索された第２の解候補から離れた（乖離度の大きい）解候補群には、同等の最が存在しないことを確認することができる。

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

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

更に他の態様に係る発明は、組合せ最適化問題における最適解を探索する最適解探索装置であって、組合せ最適化問題の解空間に属する解のうち少なくとも１つの解を第１の解候補として取得する解候補取得部と、第１の解候補からの乖離度が第１の範囲以下に収まる解候補群を二分決定グラフとして列挙索引化する列挙索引化部であって、二分決定グラフは、組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、組合せ可能パターンを縮約して列挙索引化するデータ構造を有する、列挙索引化部と、列挙索引化された解候補群から解候補群の一部又は全部を第２の解候補として一様抽出する解候補抽出部と、第１の解候補及び第２の解候補にそれぞれ評価値を付与する評価値付与部と、第１の解候補の評価値及び第２の解候補の評価値の１つ以上の評価値に基づいて第１の最適解の探索の終了の是非を判断する探索終了判断部と、第１の最適解の探索が終了していないと判断された場合は、第２の解候補の中から選択された１つ以上の解候補であって、第１の解候補と異なる解候補を第１の解候補として更新し、解候補抽出部、評価値付与部及び探索終了判断部による処理を繰り返し実行させる制御部と、第１の最適解の探索が終了したと判断された場合は、終了と判断された評価値が付与された第１の解候補を、第１の最適解として出力する出力部と、を備える。

本発明の更に他の態様に係る最適解探索装置において、１つ以上の解候補の制約条件を受け付ける制約条件受付部を備え、列挙索引化部は、第１の解候補からの乖離度が第１の範囲以下に収まり、かつ制約条件を満たす解候補群を二分決定グラフとして列挙索引化することが好ましい。

本発明によれば、広い解空間であっても効率的に解候補群を列挙索引化して最適解を検索することができ、解候補群が評価可能な個数を超えたとしても一様抽出することで均等な探索（精度の高い探索）が可能である。

ＲＮＡ発現行列データを示す図表及び遺伝子制御ネットワークを示す図 最適解の探索の手順を示す概念図 深度が浅い探索アルゴリズムでは、局所解が選択されてしまう課題を説明するために用いた図 本発明に係る最適解探索装置のハードウェア構成を示すブロック図 図４に示した最適解探索装置の機能を示す機能ブロック 要素「Ａ，Ｂ，Ｃ，Ｄ」からなる全体集合「Ｇ_∧」を示す図 ４つの部分集合（Ｇ_（１），Ｇ_（２），Ｇ_（３），Ｇ_（４））を示す図 ３つの部分集合（Ｇ_（１），Ｇ_（３），Ｇ_（４））の選択に対応する「パターン」を示す図 部分集合の選択に対応する「パターン」の３つの例と、各「パターン」が条件を満たしているか否かの判定結果とを示す図 図６から図９に示した集合分割問題における全てのパターン及び判定結果を、「２分グラフ」で網羅的に表現した図 フロンティア法の「枝刈り」により組合せ可能パターンが縮約される様子を示す図 フロンティア法の「節共有」により組合せ可能パターンが縮約される様子を示す図 フロンティア法の「枝刈り」及び「節共有」により組合せ可能パターンを縮約した結果を示す図 ノード＝｛Ａ，Ｂ，Ｃ｝で、エッジの数が３本のグラフ集合をＺＤＤ表現した図及び特定のパス（Ａ⇔Ｂ⇒Ｃ）に対応するグラフを示す図 ＺＤＤ表現したグラフ集合の全体のグラフ数（採用総数）の「数え上げ」を示す図 ＺＤＤ表現したグラフ集合から任意の指定番号のグラフを取り出す方法を示す図 ある遺伝子制御ネットワーク推定に対して最良の第２の解候補の判定例を説明するために用いたグラフ 本発明に係る最適解探索方法の実施形態を示すフローチャート 図１８に示した解候補群の列挙索引化するステップＳ１４の代わりに適用される、制約条件Ｃを考慮して解候補群の列挙索引化するステップＳ１４_1を示すフローチャート 図１８に示したステップＳ２０における第２の解候補GN_bestの第１の判定方法を示すフローチャート 図１８に示した最適解探索方法の変形例を示すフローチャート 図１８に示したステップＳ２０における第２の解候補GN_bestの第２の判定方法を示すフローチャート

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

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

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

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

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

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

本問題の解を与える手法を、図２を参照しながら説明する。

（１）組合せ最適化問題の解空間（探索空間）に属する解のうち少なくとも１つの解として適当な初期解（第１の解候補）G_0を与える。例えば、ネットワーク推定問題では空グラフG_0={}を初期解とする。また、現在解（第１の解候補）をG_xとする。最初の第１の解候補G_xは、初期解G_0である（G_x=G_0）。

（２）第１の解候補G_xからの乖離度Ｎが（第１の範囲N_1）以下に収まる解候補群GN_set={GN_1,GN_2,…GN_n}を列挙する。ここで、乖離度とは、例えば集合間のハミング距離等で定義できる。乖離度Ｎは、１以上二分決定グラフとして列挙索引化が可能な最大の乖離度以下の適当な一定値で、最も単純には乖離度Ｎ＝１としてもよい。あるいは、繰り返し毎に乖離度Ｎの大きさ変化させてもよい。例えば、乖離度Ｎとして、最初に大きな値を与え、繰り返し毎に徐々に小さい値に更新してもよい。

（３）解候補群に各々評価値Ｓ（「スコア」とも言う）を与え、第１の解候補G_x及び解候補群から最良の解候補（第２の解候補）GN_bestを選ぶ。通常、解候補群すべてに評価値Ｓを与えるが、大きな乖離度Ｎに対しては全数評価が難しい場合もある。そのような場合は、評価可能な個数Ｔ個を取り出し、Ｔ個の中から最良の解候補を選ぶ。

（４）何らかの終了条件を満たした場合は、最適解（第１の最適解）G_1=GN_bestとしてその最適解を出力する。そうでない場合は、第１の解候補G_xを乖離度Ｎ以内の最良の解候補（第２の解候補）GN_bestに更新して（G_x=GN_bestとして）上記の（１）に戻る。

尚、図２において、N_0、N_1、N_2、N_3は、それぞれ初期解G_0、第１の解候補G_xからの乖離度がＮ以内の探索範囲を示している。

上記手順に則ったアルゴリズムとして、グリーディ・ヒルクライミング法が知られている。このようなヒューリスティックな探索では探索された解が局所解に陥りやすく、特に乖離度Ｎに小さい値を与えると、局所解に陥りやすいことが知られている。

図３に示すように、乖離度が小さい（深度が浅い）探索アルゴリズムの場合、深さ方向の探索範囲が限られ、途中で探索が打ち切られることになり、近視眼的に局所解を選んでしまうという問題がある。このような局所解に陥る現象は、「水平線効果」と呼ばれている。尚、図３において、ノード間の１本の矢印は、乖離度１を示す。

これに対し、水平線の先まで見通せるような、視野の広い探索方法を構築するには、ある程度以上、大きい乖離度Ｎ以内に収まる解候補群から最適解を探索する必要がある。

しかしながら、乖離度Ｎに大きな値を与えると、局所解に陥りにくくなるものの、解候補の個数が過大になり、最適解の探索が困難になる。

例えば、ネットワーク探索問題で循環グラフを許さないベイジアンネットワーク等では、ハミング距離で乖離度を定義すると、エッジの向きを反転させるのに少なくともＮ＝２が必要である（※<(A,B),(B,A)>=<1,0>→<(A,B),(B,A)>=<0,1>）。したがって、Ｋ本のエッジが反転したグラフを探索範囲に収めたい場合、Ｎ≧２Ｋである必要が生じ、Ｋ=５ならＮ=１０等が必要となる。

第１の解候補からの乖離度が大きな解候補を、当然に明らかに不適な解（不能解）を除いて列挙（前述の手順（２）に示した解候補群GN_setを列挙）することは難しく、しばしば乖離度Ｎを十分小さい値とせざるを得ない場合がある。

例えば、以下に示す表１は、ノード数２６の適当なネットワークに対する、乖離度Ｎ（２〜１０）に対応する解候補の個数を示す表である。

表１に示すように乖離度Ｎが大きくなるにつれ、指数的に解候補の個数が増大することが分かる。これはネットワーク規模が大きくなればさらに顕著になる。

そこで、本発明は二分探索木を用いる。ここでいう二分探索木は、例えばＺＤＤ（Zero-suppressed binary Decision Diagram）であり、解を一様抽出する手段でもある。これによって極めて大きな個数の組み合わせも効率的に数えられることができる。実際、表１もＺＤＤを用いて個数をカウントしたものである。

これは、「Ｉ．組合せの一部によって、残りを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減する手段」と、「II．組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの識別工程を共有することで、識別すべきパターンを削減する手段」とを両方とも有している。

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

このうち、「Ｉ．不適かどうか」の判定には制約条件Ｃを用いることができる。

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

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

大きな乖離度Ｎに対して解を生成する方法としては、例えばランダム生成も考えられる。

また、従来のヒューリスティック探索では、例えば、乖離度Ｎ＝１等の小さな乖離度Ｎに対して、スコアＳの罰則付き最大尤度に乱数等を足して、確率的に必ずしも最適な方向に進まない手法なども提案されてきた。あるいはスコアＳに基づく確率によって次の解候補を選ぶ手法もあった。これらは、解空間のより広い範囲を探索して、局所解に陥りにくくする効果を狙ったものである。しかし、結局はＮ＝１の範囲内の解候補自体のスコアの影響が大きく、必ずしも期待通りの効果が得られるとは限らない。本発明に従うと、乖離度Ｎの実際に広い範囲内の解候補を効率的に探索できるので、より本質的な解決が期待できる。

さらに、乖離度Ｎの範囲内の解の個数が非常に大きくなり、それらを評価し切れない場合も想定される。この場合でも、本発明で提案するＺＤＤ等の二分探索木を用いることで、評価可能な個数の解を一様抽出できるので、確率的に均等に広大な解空間を探索できることになる。もちろん、本発明による大きな乖離度Ｎに対して前述の工夫を組み合わせても構わない。

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

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

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

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

［実施形態］
図５は、図４に示した最適解探索装置１０のＣＰＵ１２の機能を示す機能ブロック図である。

ＣＰＵ１２は、ハードディスク装置２０に格納された最適解探索プログラムを実行することより各種の処理部として機能し、解候補取得部１００、列挙索引化部１０２、解候補抽出部１０４、評価値付与部１０６、探索終了判断部１０８及び制御部１１０としての機能を有する。

解候補取得部１００は、組合せ最適化問題の解空間（探索空間）に属する解のうち少なくとも１つの解として適当な初期解（第１の解候補）G_0を取得する。例えば、ネットワーク推定問題では空グラフG_0={}を初期解とすることができる。また、現在解を第１の解候補G_xとして取得する。

尚、第１の解候補G_xの取得の詳細については後述する。また、最初の第１の解候補G_xは、初期解G_0を使用する（G_x=G_0）。

列挙索引化部１０２は、第１の解候補G_x（最初の第１の解候補G_xは、初期解G_0）から乖離度Ｎ以下に収まる解候補群GN_setを、二分決定グラフとして列挙索引化する部分であり、本例では、ＺＤＤを用いたパス列挙索引化アルゴリズムにより解候補群GN_setを列挙索引化する。尚、解候補群GN_setは、第１の解候補（初期解G_0又は現在の第１の解候補G_x）の探索範囲（図２のＮ_0、Ｎ_1、Ｎ_2、Ｎ_3等で示した各探索範囲）以内の解である。

制約条件Ｃの下でＺＤＤを構築することで、解候補群GN_setの総数と、集合｛Ｇ｝の任意の要素を一様抽出できるようになる。遺伝子制御ネットワークでは、遺伝子をノード、制御関係をエッジとするグラフを考える。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

図１５は、全体の解候補数（採用総数）の「数え上げ」を示す図である。

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

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

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

図１６は、任意の指定番号の解候補を取り出す方法を示す図である。

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

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

このように採用総数の「数え上げ」後、解候補取得部１００は、採用総数内の任意の番号を指定すると、指定した番号により一意に特定される解候補を取り出すことができる。これにより採用総数以下の乱数を発生させることで、解空間（初期解G_0又は現在の第１の解候補G_x）から乖離度がＮ以内の探索範囲上の解候補を「一様抽出」することができる。この解空間上の解の「一様抽出」は、ＺＤＤの重要な性質の一つである。

図５に戻って、解候補抽出部１０４は、列挙索引化部１０２と協働して、組合せ最適化問題（遺伝子制御ネットワーク）の解空間上の解（第２の解候補Ｇ）を一様抽出する部分であり、本例では、第１の解候補G_xから乖離度Ｎ（第１の範囲）以下に収まる、列挙索引化部１０２により列挙索引化された解候補群GN_set={GN_1,GN_2,…GN_n}から解候補群GN_setの一部又は全部を第２の解候補Ｇとして一様抽出する。

大きな乖離度Ｎに対して第２の解候補Ｇを一様抽出する手段としては、例えばランダム生成が考えられる。即ち、解候補をランダムに生成（例：グラフのエッジのありなしを乱数で決定）し、Ｇの制約を満たしていなければ再生成を繰り返す。しかし、ランダム生成では、既に生成済の解候補と同一の解候補が生成されてしまう可能性がある。したがって、本手法では大きな乖離度Ｎを与えると、考えられる解候補が非常に多くなり、それを網羅するか、もしくは、そこから抜き取る解の個数を大きくしないと十分な探索ができない場合がある。このような場合、ランダム生成では、抜き取り個数が大きくなれば、それだけ同一解候補を生成してしまう確率が大きくなり、したがって、解を列挙索引化できるＺＤＤ導入の効果が大きい。

さらに、特にネットワーク探索問題で循環グラフを禁止したい場合等、単純なランダム生成において制約をかけるのは難しい場合がある。特に判定・再生成を繰り返す手法で、ランダム生成がカバーする解空間に対して制約下の解空間が小さい場合、例えば解空間のサイズを１：Ｋとすれば、１個の解を生成するのに平均してＫ個のランダム生成が必要になるので、効率が悪い。乱数的な生成法ではＫが大きくなることが想定される。したがって、制約条件を課した中で解探索する場合、ＺＤＤ導入はさらに大きな効果を期待できる。

また、前述のような制約条件は、もちろんどの組み合わせ問題にも設定し得るが、特にネットワーク探索問題において、循環グラフの禁止の他、全域木（森や林）などに限る、次数に注目したスケールフリー性を担保させるなど、様々な制約が考えられるので、本発明が特に好適となる。

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

探索終了判断部１０８は、解候補取得部１００により取得した現在の第１の解候補G_xの評価値S(G_x)と、第２の解候補Ｇの中で選択された最良の第２の解候補Ｇである第２の解候補GN_bestの評価値S(GN_best)とに基づいて、最適解（第１の最適解）の探索終了を判断する。ここで、探索終了の判断としては、例えば以下の方法が考えられる。
(ア) G_x=GN_bestになる（解が更新されなくなる）
(イ) S(GN_best)-S(G_x)＜△Sになる（評価値の改善度合いが小さくなる）
(ウ) 繰り返し回数が一定回数に達する
制御部１１０は、解候補取得部１００、列挙索引化部１０２、解候補抽出部１０４、評価値付与部１０６、及び探索終了判断部１０８の各処理部を統括制御する部分であり、探索終了判断部１０８が探索終了と判断しない場合は、現在の第１の解候補G_xを、第２の解候補GN_bestに更新する（G_x=GN_best）。解候補取得部１００は、更新された新たな第１の解候補G_x、及びその評価値S(G_x)（=S(GN_best))を取得する。

また、制御部１１０は、更新された新たな第１の解候補G_xに基づいて、列挙索引化部１０２、解候補抽出部１０４、評価値付与部１０６、及び探索終了判断部１０８の各処理を繰り返し実行させる。

一方、制御部１１０は、探索終了判断部１０８が探索終了と判断した場合は、現在の第１の解候補G_xを、最適解（第１の最適解）G_1（=GN_best）として出力する。

［第２の解候補GN_bestの判定］
一様抽出された第２の解候補Ｇの中で、最大の評価値S(GN_best)が付与された最良の第２の解候補Ｇを、第２の解候補GN_bestとして抽出することができるが、第２の解候補GN_bestが、現在の第１の解候補G_xから乖離度Ｎ以下に収まる、二分決定グラフとして列挙索引化された解候補群GN_setの中で、ＲＮＡ発現行列データＤに対して真に最も適合しているかどうかを判断できないという問題がある。そのため、多額のコストを要する介入実験に踏み込むためには、例えば第１の解候補G_xをバイオロジストが精査して妥当性を判断するなどの属人的で不確実な工程を要していた。

本発明の好ましい実施形態では、最適と思われる第１の解候補G_xの最適性を、容易に見極めることができるようにする。

最適と思われる第１の解候補G_xの最適性を見極めるために、第２の解候補Ｇの抽出個数を超える個数を想定した場合の評価値のうちの最大評価値（第１の最大評価値）Ｚを推定する。

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

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

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

一様抽出した第２の解候補Ｇに対応する評価値Ｓと推定した信頼区間付の第１の最大評価値Ｚとを比較し、比較結果に基づいて第２の解候補Ｇの評価値Ｓが、第１の最大評価値Ｚの信頼区間内に入るか否かを判定する。第１の最大評価値Ｚの信頼区間内に収まっていれば、その第２の解候補Ｇは、解空間全体の最良の第２の解候補GN_bestの一つ（「第２の解候補GN_bestは十分」）であることが分かる。

仮に、Ｚ≫S(G_x)であれば、両者の評価値の差を解空間上の距離に変換し、現在推定している第２の解候補Ｇが真の第２の解候補GN_best（第１の最大評価値Ｚに対応する第２の解候補Ｇ）からどのくらい離れているかを推定してもよい。

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

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

そこで、第２の解候補GN_bestの最適十分性判定に成功した場合、第１の解候補G_xからの乖離度Ｎが、第１の範囲（乖離度N1）以内であって、第２の範囲（乖離度N2)以外の解空間（部分空間）上の第２の解候補Ｇの評価値のうちの最大評価値（第２の最大評価値）Ｗを推定する。

＜第２の最大評価値Ｗの推定＞
第２の最大評価値Ｗの推定は、第１の解候補G_xからの乖離度Ｎが、第１の範囲（乖離度N1）以内であって、第２の範囲（乖離度N2)以外の部分空間上の第２の解候補Ｇにそれぞれ対応する評価値Ｓを付与し、付与した評価値Ｓに基づいて、第２の解候補Ｇの抽出個数を超える個数を想定した場合の評価値のうちの最大評価値（第２の最大評価値Ｗ）を推定する。

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

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

そして、最適十分性判定に成功した第２の解候補GN_bestに対応する評価値S(GN_best)と、推定した信頼区間付の第２の最大評価値Ｗとを比較し、その比較結果に基づいて第２の解候補GN_bestの評価値S(GN_best)が、第２の最大評価値Ｗ（信頼区間付の第２の最大評価値Ｗの範囲）を超えているか否かを判定する。

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

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

尚、第１の解候補G_xから離間させる第２の範囲（乖離度N2）は、事前に設定する必要があるが、これはグラフの特性から設定してもよいし、経験的に設定してもよいし、例えばS(GN_best)≫Ｗになるところまで徐々に離れる乖離度N2を大きくしていってもよい。例えば、十分大きな乖離度N2の設定値から２分検索などによって効率的に適切な乖離度N2を繰り返し探索しても構わない。

言うまでもないことだが、乖離度N2がゼロであれば、Ｗ≒Ｚとなり、最適十分性判定と変わらない結果になろうし、ゼロ以外の最短の乖離度N2の場合は、第２の解候補GN_bestのみしか排除しないため、ある程度の乖離度N2を設定しないと、結果には大きな違いは出にくいと想定される。

これにより、第２の最大評価値Ｗを超えていれば、乖離度N2は、探索空間全域での唯一の第２の解候補GN_best（「第２の解候補GN_bestは必要十分」）であることが分かる。即ち、第２の解候補GN_best以外に、最良の解候補Ｇが存在しないと判断できるようになる。

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

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

図１７は、ある遺伝子制御ネットワーク推定に対して上記の第２の解候補GN_bestの判定例を説明するために用いたグラフである。

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

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

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

また、実線はヒューリスティック探索で得られた第２の解候補に対応する実測値を示すが、第２の解候補に対する第１の最大評価値Ｚ及び第２の最大評価値Ｗが正しく推定できていることを示す。

尚、第１の最大評価値Ｚの推定、第２の最大評価値Ｗの推定、及びこれらの第１の最大評価値Ｚ及び第２の最大評価値Ｗに基づく第２の解候補GN_bestの判定は、各処理部を統括制御する制御部１１０が行ってもよいし、第１の最大評価値Ｚの推定、第２の最大評価値Ｗの推定、及び第２の解候補GN_bestの判定を、それぞれ個別の処理部にて行うようにしてもよい。

上記のように推定した第１の最大評価値Ｚ及び第２の最大評価値Ｗを用いることで、第２の解候補GN_bestが、現在の第１の解候補G_xからの乖離度がＮ以下に収まる解候補群GN_setの中で、ＲＮＡ発現行列データＤに対して真に最も適合しているとの判定を行うことができ、この場合、構築するＺＤＤを、最良の第２の解候補GN_bestの探索と、第１の最大評価値Ｚ及び第２の最大評価値Ｗの推定の両方に共有できるため、特に効率がよい。

以上から、ＲＮＡ発現行列データＤを巧く説明する遺伝子制御ネットワークＧを効率的に探索することができるようになる。これによって、Ｇを巧く推定することで、患者の層別化や薬剤メカニズムの推定、それらの介入実験の意思決定がやりやすくなる。

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

図１８において、図５に示した解候補取得部１００は、まず、初期解G_0(=G_x)を取得する（ステップＳ１０、第１のステップ）。評価値付与部１０６は、取得した第１の解候補G_x（この場合、初期解G_0）に評価値S(G_x)を付与する（ステップＳ１２、第２のステップ）。

列挙索引化部１０２は、第１の解候補G_x（最初の第１の解候補G_xは、初期解G_0）からの乖離度Ｎが第１の範囲以下に収まる解候補群GN_setを、二分決定グラフとして列挙索引化する（ステップＳ１４、第３のステップ）。本例では、前述したようにＺＤＤを用いたパス列挙索引化アルゴリズムにより解候補群GN_setを列挙索引化する。

また、第１の解候補G_xから乖離度Ｎ以下に収まる解候補群GN_setを、ＺＤＤを用いて列挙索引化する場合、組合せ最適化問題に応じて何らかの制約条件Ｃを課すことが好ましい。

制約条件Ｃはモデル上、もしくは事前知識によって課すことができる。例えばベイジアンネットワークモデルでは循環グラフを表現できないため、遺伝子制御ネットワークとして表されるグラフＧは循環グラフであってはならない。また、事前知識によってスケールフリーネットワーク性が期待されていれば、そのような制約条件Ｃを設けることも考えられる。また、組合せ最適化問題によっては制約条件Ｃを特に考慮しなくてもよい場合もある。

図１９は、図１８のステップＳ１４の代わりに適用される、制約条件Ｃを考慮した場合のステップＳ１４_1に関して示している。

図１９において、最適解探索装置１０は、１つ以上の解候補の制約条件Ｃを受け付ける制約条件受付部（キーボード３０、マウス３２等）を備え、列挙索引化部１０２は、制約条件受付部から制約条件Ｃを受け付ける（ステップＳ５０、第７のステップ）。

列挙索引化部１０２は、第１の解候補G_xからの乖離度Ｎが第１の範囲以下に収まり、かつ制約条件Ｃを満たす解候補群GN_setを、ＺＤＤを用いて二分決定グラフとして列挙索引化する（ステップＳ５２）。

図１８に戻って、解候補抽出部１０４は、列挙索引化された解候補群GN_setから解候補群GN_setの一部又は全部を第２の解候補Ｇとして一様抽出する（ステップＳ１６、第４のステップ）。

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

次に、探索終了判断部１０８は、第２の解候補Ｇに付与された評価値Ｓに基づいて第２の解候補Ｇの中で最良の解候補Ｇである第２の解候補GN_bestを選択する（ステップＳ２０）。

探索終了判断部１０８は、現在の第１の解候補G_xの評価値S(G_x)、及び第２の解候補Ｇに付与された評価値Ｓのうちの１つ以上の評価値に基づいて、最適解（第１の最適解）の探索の終了の是非を判断する（ステップＳ２２、第６のステップ）。具体的には、現在の第１の解候補G_xとステップＳ２０で選択された第２の解候補GN_bestとが一致する場合（例えば、現在の第１の解候補G_xの評価値S(G_x)と、第２の解候補GN_bestに付与された評価値S(GN_best)とが一致する場合）、組合せ最適化問題の解空間全域における最適解の探索が終了したと判断する。

ステップＳ２２において、最適解の探索が終了していないと判断された場合（「No」の場合）、制御部１１０は、第２の解候補Ｇの中から選択された１つ以上の解候補であって、現在の第１の解候補G_xと異なる解候補を、次回の第１の解候補G_xとして更新し（ステップＳ２４）、ステップＳ１４に遷移させる。例えば、次回の第１の解候補G_xは、第２の解候補GN_bestに置き換えることができ、また、この場合の次回の第１の解候補G_xの評価値S(G_x)は、第２の解候補GN_bestに付与された評価値S(GN_best)とすることができる。

そして、ステップＳ２２において、最適解の探索が終了と判断されるまで、ステップＳ１４からステップＳ２４の処理が繰り返され、その結果、図２で説明したように第１の解候補G_xは、次々と評価値の高い第１の解候補G_xへと更新される。勿論、第１の解候補G_xからの乖離度がＮ以内の探索範囲（N_0、N_1、N_2、…）も次々と更新されることは言うまでもない。尚、乖離度Ｎは、一定値でもよいし、第１の解候補G_xが更新される毎に変化する値でもよい。例えば、乖離度Ｎとして最初に大きな値を与え、第１の解候補G_xが更新される毎に乖離度Ｎを徐々に小さい値に更新することが考えられる。

一方、ステップＳ２２において、最適解の探索が終了と判断された場合（「Yes」の場合）、出力部として機能する制御部１１０は、現在の第１の解候補G_xを、最適解（第１の最適解）G_1（=GN_best）として出力する（ステップＳ２６）。

尚、ステップＳ２２における最適解の探索終了の判断は、上記の実施形態に限らず、例えば、第２の解候補GN_bestの評価値S(GN_best)と現在の第１の解候補G_xの評価値S(G_x)との差が、誤差に相当する閾値△S未満になる場合（S(GN_best)-S(G_x)＜△S）とすることできる。この場合、第１の解候補G_xの評価値S(G_x)の改善度合いが小さくなり、ここで探索を打ち切っても所望の最適解が得られるからである。

また、ステップＳ１４からステップＳ２４の処理の繰り返し回数（第１の解候補G_xが更新される回数）が一定回数に達する場合を、最適解の探索終了と判断してもよい。一定回数としては、解空間全体の大きさ及び第１の解候補G_xからの乖離度Ｎの大きさにもよるが、第１の解候補G_xが収束したと見なせる回数であることが好ましい。

＜第２の解候補GN_bestの選択方法＞
次に、図１８に示したステップＳ２０での第２の解候補GN_bestの選択方法（判定方法）について説明する。

〈第２の解候補GN_bestの第１の判定方法〉
図２０は、図１８に示したステップＳ２０における第２の解候補GN_bestの第１の判定方法を示すフローチャートである。

図２０において、ステップＳ１８（図１８）にて評価値Ｓが付与された、一様抽出された第２の解候補Ｇの評価値Ｓに基づいて、第２の解候補Ｇの個数を超える個数の解Ｇを想定した場合の最大の評価値Ｓを、第１の最大評価値Ｚとして推定する（ステップＳ３０、第８のステップ）。

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

続いて、第２の解候補Ｇの評価値Ｓと推定した信頼区間付の第１の最大評価値Ｚとを比較し（ステップＳ３２）、その比較結果に基づいて第２の解候補Ｇの評価値Ｓが、信頼区間付の第１の最大評価値Ｚの範囲に入るか否かを判定する（ステップＳ３４、第９のステップ）。即ち、ある第２の解候補Ｇの評価値Ｓが、信頼区間付の第１の最大評価値の信頼区間内に入っている場合には、その第２の解候補Ｇは、ＺＤＤにより二分決定グラフとして列挙索引化された解候補群GN_setの中の最良の解（第２の解候補GN_best）としての十分条件を満たしている（最適十分性あり）と判定する。

最適十分性があると判定された最良の第２の解候補Ｇは、第２の解候補GN_bestとして出力される（ステップＳ３６）。尚、第２の解候補GN_bestに付与された評価値Ｓは、評価値S(GN_best）とされる。

このようにして真の最適値としての最適十分性を有する第２の解候補GN_bestを選択（判定）することができる。また、構築するＺＤＤを、第２の解候補GN_bestの探索と、第１の最大評価値Ｚの推定の両方に共有できるため、効率がよい。

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

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

図２１は、図１８に示した最適解探索方法の変形例を示すフローチャートである。

ヒューリスティックな探索法として、演算コストは小さい（探索時間は短い）が解の精度が低い第１の探索法と、第１の探索法よりも演算コストは大きい（探索時間は長い）が解の精度が高い第２の探索法とを準備しておく。

図２１において、最初に第１の探索法を適用する（ステップＳ６０）。図１８に示した最適解探索方法による最適解の探索を実施するステップＳ６２（特に、ステップＳ１４に相当するステップ）では、第１の探索法により探索された解を解候補群GN_setとして列挙探索化し、最適解の探索を実施する。

続いて、最適解が探索されて最適解の探索が終了したか否かを判別する（ステップＳ６４）。

最適解が探索されずに最適解の探索が終了した場合（「No」の場合）、第１の探索法を第２の探索法に切り替え（ステップＳ６６）、ステップＳ６２に遷移させる。これにより、ステップＳ６２では、第２の探索法により探索された解を解候補群GN_setとして列挙探索化し、最適解の探索を実施することになる。

最適解が探索されて最適解の探索が終了した場合（「Yes」の場合）、本探索が終了する。

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

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

尚、第１の探索法を使用して最適解の探索が終了した場合には、探索法の切り替えは行われない。

〈第２の解候補GN_bestの第２の判定方法〉
図２２は、図１８に示したステップＳ２０における第２の解候補GN_bestの第２の判定方法を示すフローチャートであり、特に図１８に示したステップＳ２２において、最適解の探索が終了したと判断された後（即ち、最適十分性判定に成功した後）に行われる最適必要性判定の処理に関して示している。

図２２において、現在の第１の解候補G_xからの乖離度が、第１の範囲(乖離度N1）以下に収まる解候補群に含まれる第２の解候補群であって、第１の最大評価値Ｚの信頼区間内に評価値Ｓが入ると判定された第２の解候補Ｇ（即ち、最適十分性判定に成功した第２の解候補GN_best）からの乖離度が、第２の範囲（乖離度N2)外の第２の解候補群GN_setを二分決定グラフとして列挙索引化する（ステップＳ７０、第１０のステップ）。尚、乖離度N2は、乖離度N1よりも小さい値（N2＜N1）である。また、第２の解候補群GN_setは、ＺＤＤにより二分決定グラフとして列挙索引化されたものである。

続いて、列挙索引化された第２の解候補群GN_setから第２の解候補群GN_setの一部又は全部を第３の解候補Ｇとして一様抽出する（ステップＳ７２、第１１のステップ）。

一様抽出した第３の解候補Ｇに評価値Ｓを付与する（ステップＳ７４、第１２のステップ）。第３の解候補Ｇに対する評価値Ｓの付与は、図１８に示したステップＳ１８と同様に行うことができる。

次に、第３の解候補Ｇに付与された評価値Ｓに基づいて、第３の解候補の個数を超える個数の解を想定した場合の最大の評価値Ｓを、第２の最大評価値Ｗとして推定する（ステップＳ７６、第１３のステップ）。第２の最大評価値Ｗの推定は、図２０に示したステップＳ３０と同様に行うことができる。

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

続いて、第１の最大評価値Ｚの信頼区間内に入った第２の解候補GN_bestの評価値S(GN_best)と第２の最大評価値Ｗとを比較し（ステップＳ７８）、その比較結果に基づいて第２の解候補GN_bestの評価値S(GN_best）が、第２の最大評価値Ｗを超えているか否かを判定する（ステップＳ８０、第１４のステップ）。その比較結果に基づいて第２の解候補GN_bestの評価値S(GN_best）が、信頼区間付の第２の最大評価値Ｗの範囲成功した第２の解候補GN_bestに対応する評価値S(GN_best)と、推定した信頼区間付の第２の最大評価値Ｗとを比較し、その比較結果に基づいて第２の解候補GN_bestの評価値S(GN_best)が、第２の最大評価値Ｗ（信頼区間付の第２の最大評価値Ｗの範囲）を超えているか否かを判定する。

第２の解候補GN_bestの評価値S(GN_best)が、第２の最大評価値Ｗを超えている場合には、第２の解候補GN_bestと同等解は他には存在せず、第２の解候補GN_bestは唯一の最適解としての必要条件を満たしている（最適必要性あり）と判定することができる。

〈第２の解候補GN_bestの第２の判定方法の変形例〉
第２の解候補GN_bestの第２の判定方法の変形例は、図２２に示した第２の解候補GN_bestの第２の判定方法において、第２の解候補GN_bestの最適必要性判定に失敗した場合の処理を含む。

即ち、最適必要性判定に失敗した場合、図２２に示したステップＳ７０において、第２の範囲（乖離度N2)に代えて第２の範囲を拡大した第３の範囲（乖離度N3：N3＞N2）を適用して、図２２に示したステップＳ７０からステップＳ８０の処理を行う。

更に最適必要性判定に失敗した場合には、必要性判定に成功するまで第３の範囲を徐々に拡大して、図２２に示したステップＳ７０からステップＳ８０の処理を複数回繰り返してもよい。

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

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

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

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

また、本発明は、プロセッサを有する最適解探索装置であって、プロセッサが、組合せ最適化問題の解空間に属する解のうち少なくとも１つの解を第１の解候補として取得し、第１の解候補に評価値を付与し、第１の解候補からの乖離度が第１の範囲以下に収まる解候補群を、ＺＤＤ又はこれに類似するデータ構造を有する二分決定グラフとして列挙索引化し、列挙索引化された解候補群から解候補群の一部又は全部を第２の解候補として一様抽出し、抽出した第２の解候補に評価値を付与し、第１の解候補の評価値及び第２の解候補の評価値の１つ以上の評価値に基づいて第１の最適解の探索の終了の是非を判断する。第１の最適解の探索が終了していないと判断された場合は、第２の解候補の中から選択された１つ以上の解候補であって、第１の解候補と異なる解候補を第１の解候補として更新して、第１の最適解の探索が終了したと判断されるまで上記の処理を繰り返し、第１の最適解の探索が終了したと判断された場合は、終了と判断された評価値が付与された第１の解候補を、第１の最適解として出力する最適解探索装置を含む。

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

１０ 最適解探索装置
１２ ＣＰＵ
１４ 主メモリ
１６ グラフィックボード
１８ 通信インターフェース
２０ ハードディスク装置
２２ 光学ディスクドライブ
２４ キーボードコントローラ
２６ マウスコントローラ
２８ モニタ装置
３０ キーボード
３２ マウス
４０ データベース
５０ ネットワーク
１００ 解候補取得部
１０２ 列挙索引化部
１０４ 解候補抽出部
１０６ 評価値付与部
１０８ 探索終了判断部
１１０ 制御部
Ｚ 第１の最大評価値
Ｗ 第２の最大評価値

Claims (16)

1. 組合せ最適化問題における最適解をコンピュータにより探索する最適解探索方法であって、
   前記組合せ最適化問題の解空間に属する解のうち少なくとも１つの解を第１の解候補として取得する第１のステップと、
   前記第１の解候補に評価値を付与する第２のステップと、
   前記第１の解候補からの乖離度が第１の範囲以下に収まる解候補群を二分決定グラフとして列挙索引化する第３のステップであって、前記二分決定グラフは、前記組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び前記組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、前記組合せ可能パターンを縮約して列挙索引化するデータ構造を有する、前記第３のステップと、
   前記列挙索引化された解候補群から前記解候補群の一部又は全部を第２の解候補として一様抽出する第４のステップと、
   抽出した前記第２の解候補に評価値を付与する第５のステップと、
   前記第１の解候補の評価値及び前記第２の解候補の評価値の１つ以上の評価値に基づいて第１の最適解の探索の終了の是非を判断する第６のステップと、を含み、
   前記第１の最適解の探索が終了していないと判断された場合は、前記第２の解候補の中から選択された１つ以上の解候補であって、前記第１の解候補と異なる解候補を前記第１の解候補として更新し、前記第３のステップから前記第６のステップの処理を繰り返し、
   前記第１の最適解の探索が終了したと判断された場合は、前記終了と判断された評価値が付与された前記第１の解候補を、前記第１の最適解として出力する、最適解探索方法。
2. １つ以上の解候補の制約条件を受け付ける第７のステップを含み、
   前記第３のステップは、前記第１の解候補からの乖離度が前記第１の範囲以下に収まり、かつ前記制約条件を満たす前記解候補群を二分決定グラフとして列挙索引化する、
   請求項１に記載の最適解探索方法。
3. 前記第１の範囲の乖離度は、１以上前記二分決定グラフとして列挙索引化が可能な最大の乖離度以下である請求項１又は２に記載の最適解探索方法。
4. 前記第１の範囲の乖離度は、一定値又は前記第１の解候補が更新される毎に変化する値である請求項３に記載の最適解探索方法。
5. 前記第６のステップは、前記第１の解候補の評価値が全ての前記第２の解候補の評価値以上の場合、前記第２の解候補と前記第１の解候補との差が規定値以下の場合、又は前記第３のステップから前記第６のステップの処理の繰り返し回数が一定回数に達する場合を、前記第１の最適解の探索の終了と判断する、
   請求項１から４のいずれか１項に記載の最適解探索方法。
6. 前記第１の解候補として更新される前記第１の解候補と異なる解候補は、前記第２の解候補の中で最大の評価値が付与された第２の解候補である、請求項１から５のいずれか１項に記載の最適解探索方法。
7. 前記第６のステップは、
   前記第５のステップで付与した前記第２の解候補の評価値に基づいて、前記第２の解候補の個数を超える個数の解を想定した場合の最大の評価値を、第１の最大評価値として推定する第８のステップと、
   前記第２の解候補の評価値が前記第１の最大評価値の信頼区間内に入るか否かを判定する第９のステップと、を含み、
   前記第１の解候補の評価値及び前記第１の最大評価値の信頼区間内に入ると判定された前記第２の解候補の評価値に基づいて前記第１の最適解の探索の終了の是非を判断する、
   請求項１から６のいずれか１項に記載の最適解探索方法。
8. 前記一様抽出される前記第２の解候補は、Ｕ，Ｖをそれぞれ自然数とすると、Ｕ×Ｖ個の解であり、
   前記第８のステップは、前記Ｕ×Ｖ個の解をＶ個のブロックに分け、前記ブロック毎にＵ個の解の評価値の区分最大値をＶ個取得し、前記Ｖ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第１の最大評価値を推定する、
   請求項７に記載の最適解探索方法。
9. 演算コストは小さいが解の精度が低い第１の探索法と、前記第１の探索法よりも演算コストは大きいが解の精度が高い第２の探索法とを有し、
   前記第３のステップは、最初に前記第１の探索法により探索された解を前記解候補群として列挙索引化し、前記第２の解候補の評価値が前記第１の最大評価値の信頼区間内に入らない場合のみ、前記第２の探索法により探索された解を前記解候補群として列挙索引化する、
   請求項７又は８に記載の最適解探索方法。
10. 前記第１の解候補からの乖離度が第１の範囲以下に収まる解候補群に含まれる第２の解候補群であって、前記第１の最大評価値の信頼区間内に評価値が入ると判定された前記第２の解候補からの乖離度が、前記第１の範囲よりも狭い第２の範囲外の前記第２の解候補群を二分決定グラフとして列挙索引化する第１０のステップであって、前記二分決定グラフは、前記組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び前記組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、前記組合せ可能パターンを縮約して列挙索引化するデータ構造を有する、前記第１０のステップと、
    前記列挙索引化された第２の解候補群から前記第２の解候補群の一部又は全部を第３の解候補として一様抽出する第１１のステップと、
    抽出した前記第３の解候補に評価値を付与する第１２のステップと、
    前記第１２のステップで付与した前記第３の解候補の評価値に基づいて、前記第３の解候補の個数を超える個数の解を想定した場合の最大の評価値を、第２の最大評価値として推定する第１３のステップと、
    前記第１の最大評価値の信頼区間内に入った前記第２の解候補の評価値が、前記第２の最大評価値を超えているか否かを判定する第１４のステップと、
    を含む請求項７から９のいずれか１項に記載の最適解探索方法。
11. 前記一様抽出される前記第３の解候補は、Ｐ，Ｑをそれぞれ自然数とすると、Ｐ×Ｑ個の解であり、
    前記第１３のステップは、前記Ｐ×Ｑ個の解をＱ個のブロックに分け、前記ブロック毎のＰ個の解の評価値の区分最大値をＱ個取得し、前記Ｑ個の区分最大値を用いて、前記区分最大値が一般極値分布に従うものとして前記第２の最大評価値を推定する、
    請求項１０に記載の最適解探索方法。
12. 前記第１４のステップにより前記第１の最大評価値の信頼区間内に入った前記第２の解候補の評価値が、前記第２の最大評価値を超えていないと判定されると、前記第２の範囲に代えて前記第２の範囲を拡大した第３の範囲を適用して前記第１０のステップから前記第１４のステップの処理を行う、
    請求項１０又は１１に記載の最適解探索方法。
13. 組合せ最適化問題は、遺伝子制御ネットワークの組合せ最適化問題である請求項１から１２のいずれか１項に記載の最適解探索方法。
14. 請求項１から１３のいずれか１項に記載の最適解探索方法をコンピュータに実行させる最適解探索プログラム。
15. 組合せ最適化問題における最適解を探索する最適解探索装置であって、
    前記組合せ最適化問題の解空間に属する解のうち少なくとも１つの解を第１の解候補として取得する解候補取得部と、
    前記第１の解候補からの乖離度が第１の範囲以下に収まる解候補群を二分決定グラフとして列挙索引化する列挙索引化部であって、前記二分決定グラフは、前記組合せ最適化問題における組合せ可能なパターンのうち、組合せの一部により残りの組合せを考慮せずとも不適となることが確定するかどうかを識別することで、識別すべきパターンを削減するステップ、及び前記組合せ可能なパターンのうち、組合せの一部だけに差分があるパターン群の共通部分を抽出し、残りの組合せを共有することで、識別すべきパターンを削減するステップのうちの少なくとも一方を用いて、前記組合せ可能パターンを縮約して列挙索引化するデータ構造を有する、前記列挙索引化部と、
    前記列挙索引化された解候補群から前記解候補群の一部又は全部を第２の解候補として一様抽出する解候補抽出部と、
    前記第１の解候補及び前記第２の解候補にそれぞれ評価値を付与する評価値付与部と、
    前記第１の解候補の評価値及び前記第２の解候補の評価値の１つ以上の評価値に基づいて第１の最適解の探索の終了の是非を判断する探索終了判断部と、
    前記第１の最適解の探索が終了していないと判断された場合は、前記第２の解候補の中から選択された１つ以上の解候補であって、前記第１の解候補と異なる解候補を前記第１の解候補として更新し、前記解候補抽出部、前記評価値付与部及び前記探索終了判断部による処理を繰り返し実行させる制御部と、
    前記第１の最適解の探索が終了したと判断された場合は、前記終了と判断された評価値が付与された前記第１の解候補を、前記第１の最適解として出力する出力部と、
    を備えた最適解探索装置。
16. １つ以上の解候補の制約条件を受け付ける制約条件受付部を備え、
    前記列挙索引化部は、前記第１の解候補からの乖離度が前記第１の範囲以下に収まり、かつ前記制約条件を満たす前記解候補群を二分決定グラフとして列挙索引化する、
    請求項１５に記載の最適解探索装置。

Images