知財求人 - 知財ポータルサイト「IP Force」
特許7219402最適化装置、最適化装置の制御方法及び最適化装置の制御プログラム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-01-31
(45)【発行日】2023-02-08
(54)【発明の名称】最適化装置、最適化装置の制御方法及び最適化装置の制御プログラム
(51)【国際特許分類】
G06N 99/00 20190101AFI20230201BHJP
【ＦＩ】
G06N99/00 180
【請求項の数】
9
(21)【出願番号】P 2019119604
(22)【出願日】2019-06-27
(65)【公開番号】P2021005282
(43)【公開日】2021-01-14
【審査請求日】2022-03-08
(73)【特許権者】
【識別番号】000005223
【氏名又は名称】富士通株式会社
(74)【代理人】
【識別番号】110002918
【氏名又は名称】弁理士法人扶桑国際特許事務所
(72)【発明者】
【氏名】宮澤 俊之
(72)【発明者】
【氏名】柴▲崎▼ 崇之
(72)【発明者】
【氏名】上村 泰紀
【審査官】加藤 優一
(56)【参考文献】
【文献】特開２０１６－０５１３５０（ＪＰ，Ａ）
【文献】特開平１１－２１２９４８（ＪＰ，Ａ）
【文献】特開平０７－１５２７７５（ＪＰ，Ａ）
【文献】株式会社富士通研究所，中分子創薬へ適用可能な組合せ最適化問題を解く技術を開発，プレスリリース ［オンライン］，日本，富士通株式会社，2018年09月18日，［検索日 2022.12.20］，インターネット: ＜ＵＲＬ：https://pr.fujitsu.com/jp/news/2018/09/18.html＞
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｎ ３／００ －９９／００
Ｇ１６Ｂ ５／００ －９９／００
Ｇ０６Ｑ １０／００ －９９／００
(57)【特許請求の範囲】
【請求項1】
組合せ最適化問題を表す第１の評価関数に含まれる複数の状態変数のそれぞれの間の相互作用の大きさを表す重み係数を保持する記憶部と、前記第１の評価関数に基づいて、前記組合せ最適化問題の局所解を算出する局所解算出部と、
前記複数の状態変数のうち前記局所解における値が所定値である第１の状態変数と、前記重み係数と、入力される選択領域情報とに基づいて、前記複数の状態変数から部分問題の対象となる状態変数群を選択する選択部と、
選択された前記状態変数群についての前記部分問題を表す第２の評価関数について、基底状態を探索する探索部と、
前記基底状態の探索によって得られた部分解に基づいて前記局所解を更新することで、全体解を生成する全体解生成部と、
を有する最適化装置。
【請求項2】
前記選択部は、前記重み係数に基づいて、前記第１の状態変数との間で相互作用がある第２の状態変数を前記状態変数群に含む候補として選択する、請求項１に記載の最適化装置。
【請求項3】
前記選択部は、前記全体解生成部が生成した前記全体解における値が前記所定値である第３の状態変数と、前記重み係数と前記選択領域情報とに基づいて、選択する前記状態変数群を更新する、請求項１または２に記載の最適化装置。
【請求項4】
前記局所解算出部は、前記全体解生成部が生成した前記全体解を初期解として、前記局所解を算出する、請求項１乃至３の何れか一項に記載の最適化装置。
【請求項5】
前記選択領域情報は、前記探索部が処理可能な演算ビット数の情報を含み、前記選択部は、前記演算ビット数と同じ数の状態変数による前記状態変数群を選択する、請求項１乃至４の何れか一項に記載の最適化装置。
【請求項6】
前記組合せ最適化問題は、複数のアミノ酸による直鎖構造の安定構造を探索する問題であり、前記複数の状態変数のそれぞれは、前記複数のアミノ酸の何れかが配置される候補となる複数の格子点の何れかに割当てられており、前記選択部は、アミノ酸が配置されている格子点に割当てられている前記第１の状態変数との間で相互作用がある前記第２の状態変数を前記状態変数群に含む候補として選択する、
請求項２に記載の最適化装置。
【請求項7】
前記選択部は、２以上のアミノ酸の配置が禁止されている同一シェルに属す格子点に割当てられた状態変数の中から、前記状態変数群を選択する、請求項６に記載の最適化装置。
【請求項8】
最適化装置が有する局所解算出部が、組合せ最適化問題を表す第１の評価関数に基づいて、前記組合せ最適化問題の局所解を算出し、前記最適化装置が有する選択部が、前記第１の評価関数に含まれる複数の状態変数のうち前記局所解における値が所定値である第１の状態変数と、記憶部に保持されている前記複数の状態変数のそれぞれの間の相互作用の大きさを表す重み係数と、入力される選択領域情報とに基づいて、前記複数の状態変数から部分問題の対象となる状態変数群を選択し、
前記最適化装置が有する探索部が、選択された前記状態変数群についての前記部分問題を表す第２の評価関数について、基底状態を探索し、
前記最適化装置が有する全体解生成部が、前記基底状態の探索によって得られた部分解に基づいて前記局所解を更新することで、全体解を生成する、
最適化装置の制御方法。
【請求項9】
最適化装置の制御プログラムにおいて、組合せ最適化問題を表す第１の評価関数に基づいて、前記組合せ最適化問題の局所解を算出し、
前記第１の評価関数に含まれる複数の状態変数のうち前記局所解における値が所定値である第１の状態変数と、記憶部に保持されている前記複数の状態変数のそれぞれの間の相互作用の大きさを表す重み係数と、入力される選択領域情報とに基づいて、前記複数の状態変数から部分問題の対象となる状態変数群を選択し、
選択された前記状態変数群についての前記部分問題を表す第２の評価関数について、基底状態を探索し、
前記基底状態の探索によって得られた部分解に基づいて前記局所解を更新することで、全体解を生成する、
処理をコンピュータに実行させる最適化装置の制御プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムに関する。
【背景技術】
【0002】
組合せ最適化問題は、現在の社会における様々な分野に存在する。たとえば、製造・流通、マーケティングなどの分野では、コストを最適化（最小化）する要素の組合せが探索される。しかし、組合せ最適化問題は、上記要素に対応する変数の数が増えるにつれて指数関数的に計算時間が増加するため、従来のノイマン型コンピュータでは解くことが困難である問題として知られている。
【0003】
このようなノイマン型コンピュータが不得意とする多変数の組合せ最適化問題を解く手法として、計算対象の組合せ最適化問題を、磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えて計算する手法がある。イジングモデルを利用して組合せ最適化問題を解くハードウェアとして、デジタル回路を用いてシミュレーテッド・アニーリングを行うことでイジング型の評価関数の値が最小となる状態変数の値の組合せを計算するものがある。また、超伝導回路を用いて量子アニーリングを行うことで同様の計算を行うものなどがある。
【0004】
一方、ハードウェア量やソフトウェアにて計算する場合のメモリサイズなどの制約により、問題規模（状態変数の数（イジングモデルのスピン数に相当））が大きい組合せ最適化問題を、小さい問題に分割して計算を行う手法がある。たとえば計算対象の組合せ最適化問題を分割した部分問題を作成し、各部分問題の解を探索し、各部分問題の解に基づいて全体の解を生成する手法があった（たとえば、特許文献１，２、非特許文献１参照）。
【先行技術文献】
【特許文献】
【0005】
【文献】特開２０１６－５１３５０号公報
【文献】特開２０１５－１９１３４０号公報
【非特許文献】
【0006】
【文献】G. Rosenberg, and et. al., “Building an iterative heuristic solver for a quantum annealer”, Computational Optimization and Applications, 65, 3, 845, 2016
【発明の概要】
【発明が解決しようとする課題】
【0007】
計算対象の組合せ最適化問題を分割して計算する場合、部分問題の解に基づいて全体の解を求める計算処理などの影響により、問題全体の解を得るまでの計算量が分割をしないで計算する場合よりも増加する問題があった。
【0008】
１つの側面では、本発明は、組合せ最適化問題を分割して計算する際の計算量の増加を抑制可能な最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムを提供することを目的とする。
【課題を解決するための手段】
【0009】
１つの実施態様では、組合せ最適化問題を表す第１の評価関数に含まれる複数の状態変数のそれぞれの間の相互作用の大きさを表す重み係数を保持する記憶部と、前記第１の評価関数に基づいて、前記組合せ最適化問題の局所解を算出する局所解算出部と、前記複数の状態変数のうち前記局所解における値が所定値である第１の状態変数と、前記重み係数と、入力される選択領域情報とに基づいて、前記複数の状態変数から部分問題の対象となる状態変数群を選択する選択部と、選択された前記状態変数群についての前記部分問題を表す第２の評価関数について、基底状態を探索する探索部と、前記基底状態の探索によって得られた部分解に基づいて前記局所解を更新することで、全体解を生成する全体解生成部と、を有する最適化装置が提供される。
【0010】
また、１つの実施態様では、最適化装置の制御方法が提供される。
また、１つの実施態様では、最適化装置の制御プログラムが提供される。
また、１つの実施態様では、最適化装置の制御プログラムが提供される。
【発明の効果】
【0011】
１つの側面では、組合せ最適化問題を分割して計算する際の計算量の増加を抑制できる。
【図面の簡単な説明】
【0012】
【図1】第１の実施の形態の最適化装置の一例を示す図である。
【図2】第１の実施の形態の最適化装置の制御方法の一例の流れを示すフローチャートである。
【図3】第２の実施の形態の最適化装置のハードウェアの一例を示す図である。
【図4】第２の実施の形態の最適化装置のハードウェアの他の例を示す図である。
【図5】第２の実施の形態の最適化装置の機能例を示すブロック図である。
【図6】組合せ最適化手法を利用したアミノ酸の安定構造の決定過程の一例を示す図である。
【図7】ダイヤモンド・エンコーディング法の例を示す図である。
【図8】各シェルの格子点に割当てられた状態変数の一例を示す図である。
【図9】各シェルにおけるアミノ酸の配置可能位置の一例を示す図である。
【図10】アミノ酸が配置されている格子点の例を示す図である。
【図11】部分問題の作成処理の一例の流れを示すフローチャートである。
【図12】保持される識別番号ｊの一例を示す図である。
【図13】部分問題の作成処理の他の例の流れを示すフローチャートである。
【図14】問題の分割の有無による計算量の比較例を示す図である。
【図15】第３の実施の形態における部分問題の作成処理の一例の流れを示すフローチャートである。
【図16】あるシェルにおいて選択候補となる状態変数の一例を示す図である。
【発明を実施するための形態】
【0013】
【0014】
第１の実施の形態の最適化装置１０は、記憶部１１、局所解算出部１２、選択部１３、探索部１４、全体解生成部１５を有する。
記憶部１１は、組合せ最適化問題を表す評価関数に含まれる複数の状態変数のそれぞれの間の相互作用の大きさを表す重み係数を保持する。記憶部１１は、たとえば、ＲＡＭ（Random Access Memory）などの揮発性の記憶装置、または、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性の記憶装置である。
記憶部１１は、組合せ最適化問題を表す評価関数に含まれる複数の状態変数のそれぞれの間の相互作用の大きさを表す重み係数を保持する。記憶部１１は、たとえば、ＲＡＭ（Random Access Memory）などの揮発性の記憶装置、または、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性の記憶装置である。
【0015】
計算対象の組合せ最適化問題としては、たとえば、複数のアミノ酸（アミノ酸残基）による直鎖構造の安定構造を探索する問題、巡回セールスマン問題などである。これらの問題は、より実用的な問題とした場合、状態変数の数が極めて大きくなることが知られており、ハードウェア量やソフトウェアにて計算する場合のメモリサイズなどの制約から、部分問題に分割して計算することが好ましい。ただ、計算対象の問題としては、これらの問題に限定されるものではない。
【0016】
重み係数を用いたイジング型の評価関数（エネルギー関数とも呼ばれる）Ｅ（ｘ）は、たとえば、以下の式（１）で定義される。
【0017】
【数1】
【0018】
右辺の１項目は、評価関数に含まれる全状態変数から選択可能な２つの状態変数の全組合せについて、漏れと重複なく、２つの状態変数の値（０または１）と重み係数との積を積算したものである。ｘｉは、識別番号がｉの状態変数、ｘｊは、識別番号がｊのビットの値を表す状態変数であり、Ｗｉｊは、ｘｉとｘｊの間の相互作用の大きさを示す重み係数である。なお、Ｗｉｉ＝０である。また、Ｗｉｊ＝Ｗｊｉであることが多い（つまり、重み係数による係数行列は対称行列である場合が多い）。
【0019】
たとえば、図１に示すように、状態変数の総数がｎ個の場合、記憶部１１は、ｎ×ｎ個の重み係数を保持している。各重み係数の値は、計算対象の組合せ最適化問題の特徴を反映したものとなっている。
【0020】
ｎ×ｎ個の重み係数は、最適化装置１０の外部から取得してもよいし、最適化装置１０が、入力された問題データに基づいて算出するようにしてもよい。
局所解算出部１２は、式（１）に示したような評価関数に基づいて、組合せ最適化問題の局所解を算出する。局所解算出部１２は、たとえば、タブーサーチ法などによって、局所解を算出する。局所解は、後述する探索部１４による探索処理のための初期解として用いられるため、最適解でなくてもよい。なお、局所解算出部１２は、後述する全体解生成部１５が生成した全体解を初期解として局所解を算出してもよい。
局所解算出部１２は、式（１）に示したような評価関数に基づいて、組合せ最適化問題の局所解を算出する。局所解算出部１２は、たとえば、タブーサーチ法などによって、局所解を算出する。局所解は、後述する探索部１４による探索処理のための初期解として用いられるため、最適解でなくてもよい。なお、局所解算出部１２は、後述する全体解生成部１５が生成した全体解を初期解として局所解を算出してもよい。
【0021】
選択部１３は、複数の状態変数のうち局所解における値が所定値である状態変数と、重み係数と、入力される選択領域情報とに基づいて、複数の状態変数から部分問題の対象となる状態変数群を選択する。
【0022】
たとえば、式（１）において、値が１である状態変数を第１の状態変数、第１の状態変数と相互作用のある状態変数を第２の状態変数とすると、第２の状態変数の値に応じて、評価関数の値が変化する可能性がある。つまり、評価関数の値が小さくなるように（より基底状態に近づくように）、第２の状態変数の値が更新される可能性がある。そこで、選択部１３は、重み係数に基づいて、第２の状態変数を状態変数群に含む候補として選択する。
【0023】
ただ、重み係数の値は、後述する評価関数の制約項を表現するものである場合、非常に大きな値になる可能性がある。そのような重み係数によって第１の状態変数との間の相互作用の大きさが示されている第２の状態変数については、値が１であると制約項が表す制約条件を満たさなくなる可能性があるため、選択部１３は、状態変数群には含めない。このような機能を実現するために、状態変数群を選択するために用いる重み係数の値の範囲などが、入力される選択領域情報に含まれる。
【0024】
また、状態変数群に含まれる状態変数の数を増やすために、たとえば、選択部１３は、第１の状態変数に対して識別番号が近い所定数の状態変数に対しても、同様に重み係数に基づいて、状態変数群に含める第２の状態変数を選択する処理を行ってもよい。上記所定数についても、たとえば、選択領域情報に含まれる。上記所定数は、たとえば、探索部１４が処理可能な演算ビット数に応じて決定される。演算ビット数についても、選択領域情報に含まれていてもよい。探索部１４による処理の効率化のために、探索部１４が処理可能な演算ビット数と同じ数の状態変数が、状態変数群に含まれることが望ましい。
【0025】
なお、選択部１３は、後述する全体解生成部１５が生成した全体解における値が所定値である状態変数と、重み係数と選択領域情報とに基づいて、選択する状態変数群を更新してもよい。
【0026】
探索部１４は、選択部１３によって選択された状態変数群についての部分問題を表す評価関数について、基底状態を探索する。探索部１４は、シミュレーテッド・アニーリング法やレプリカ交換法などのマルコフ連鎖モンテカルロ法により基底状態の探索を行ってもよいし、量子アニーリング法により基底状態の探索を行ってもよい。部分問題を表す評価関数は、式（１）に含まれる状態変数を、選択された状態変数群としたものである。
【0027】
なお、基底状態の探索によって得られる部分問題の解（部分解）である状態（選択された状態変数群の値）は、たとえば、所定回数更新された状態のうち、部分問題を表す評価関数の値が最小となる状態である。
【0028】
全体解生成部１５は、部分解に基づいて局所解を更新することで全体解を生成する。たとえば、全体解生成部１５は、部分解が示す選択された状態変数群の値を、局所解における対応する状態変数群の値と入れ替えることで全体解を生成する。
【0029】
局所解算出部１２、選択部１３、探索部１４及び全体解生成部１５は、たとえば、ＣＰＵ（Central Processing Unit）やＤＳＰ（Digital Signal Processor）などのプロセッサが実行するプログラムモジュールを用いて実装できる。なお、探索部１４は、デジタル回路を用いてシミュレーテッド・アニーリング法やレプリカ交換法などを実行するハードウェアであってもよいし、量子アニーリングを行うハードウェアであってもよい。
【0030】
以下第１の実施の形態の最適化装置１０の動作（最適化装置１０の制御方法）の例を説明する。
図２は、第１の実施の形態の最適化装置の制御方法の一例の流れを示すフローチャートである。
図２は、第１の実施の形態の最適化装置の制御方法の一例の流れを示すフローチャートである。
【0031】
最適化装置１０の図示しない入力部は、前述したような重み係数や選択領域情報を取得する（ステップＳ１）。重み係数は記憶部１１に記憶される（選択領域情報も記憶部１１に記憶されてもよい）。局所解算出部１２は、たとえば、式（１）に示したような評価関数に基づいて、組合せ最適化問題の局所解を算出する（ステップＳ２）。図１には、局所解の一例が示されている。図１の例では、ｎ個の状態変数（ｘ１～ｘｎ）のうち、ｘｉの値が１となっている。
【0032】
選択部１３は、複数の状態変数のうち局所解における値が所定値である状態変数と、重み係数と、入力される選択領域情報とに基づいて、複数の状態変数から部分問題の対象となる状態変数群を選択する（ステップＳ３）。図１には、値が１である状態変数（前述の第１の状態変数）と、重み係数と、入力される選択領域情報とに基づいて状態変数群が選択される例が示されている。たとえば、ｘｉの値が１であり、Ｗｉｊの値が０ではない（かつ、選択領域情報で表される所定の範囲内の値である）場合、ｘｊが優先的に、状態変数群の１つとして選択される。さらに、選択部１３は、ｘｉに対して識別番号が近い所定数の状態変数に対しても、同様に重み係数に基づいて、状態変数群に含める状態変数を選択する処理を行う。
【0033】
その後、探索部１４は、選択部１３によって選択された状態変数群についての部分問題を表す評価関数について、基底状態を探索することによって部分解を計算する（ステップＳ４）。そして、全体解生成部１５は、部分解に基づいて局所解を更新することで全体解を生成する（ステップＳ５）。
【0034】
その後、たとえば、全体解生成部１５は、ステップＳ３～Ｓ５の処理の反復回数（サブループ反復回数）が所定数（ｒ１）であるか否かを判定する（ステップＳ６）。サブループ反復回数＝ｒ１でないと判定された場合には、ステップＳ３からの処理が繰り返される。なお、サブループの２回目以降の処理において、選択部１３は、全体解生成部１５が生成した全体解に基づいて、部分問題の対象となる状態変数群を更新する。
【0035】
全体解生成部１５は、サブループ反復回数＝ｒ１であると判定した場合、ステップＳ２～Ｓ６の処理の反復回数（メインループ反復回数）が所定数（ｒ２）であるか否かを判定する（ステップＳ７）。メインループ反復回数＝ｒ２でないと判定された場合には、ステップＳ２からの処理が繰り返される。なお、メインループの２回目以降の処理において、局所解算出部１２は、全体解生成部１５が生成した全体解を初期解として、局所解を算出する。
【0036】
全体解生成部１５は、メインループ反復回数＝ｒ２であると判定した場合、たとえば、これまでに得られた全体解のうち、評価関数の値が最小となる全体解を、最良の全体解として出力する（ステップＳ８）。
【0037】
なお、処理の順序は、図２の例に限定されるわけではなく、適宜入れ替えてもよい。
以上のように、最適化装置１０は、組合せ最適化問題を部分問題に分割する際、局所解において値が所定値である状態変数と重み係数に基づき部分問題の対象となる状態変数群を選択することで、組合せ最適化問題の特徴を反映した適切な分割が可能になる。前述のように一般的に計算対象の組合せ最適化問題を分割して計算する場合、部分問題の解に基づいて全体の解を求める計算処理などの影響により、問題全体の解を得るまでの計算量が分割をしないで計算する場合よりも増加する。部分問題の対象となる状態変数群を、単に乱数に基づいてランダムに選択した場合、部分問題に組合せ最適化問題の特徴が反映されにくい。そのため、そのような部分問題を用いて計算を行う場合、精度のよい全体解を得るため多くの反復計算を要し、計算量がより増加する可能性がある。一方、最適化装置１０では、上記のように組合せ最適化問題の特徴を反映した適切な分割が可能になるため、精度のよい全体解を得るまでの計算量の増加を抑制できる。
以上のように、最適化装置１０は、組合せ最適化問題を部分問題に分割する際、局所解において値が所定値である状態変数と重み係数に基づき部分問題の対象となる状態変数群を選択することで、組合せ最適化問題の特徴を反映した適切な分割が可能になる。前述のように一般的に計算対象の組合せ最適化問題を分割して計算する場合、部分問題の解に基づいて全体の解を求める計算処理などの影響により、問題全体の解を得るまでの計算量が分割をしないで計算する場合よりも増加する。部分問題の対象となる状態変数群を、単に乱数に基づいてランダムに選択した場合、部分問題に組合せ最適化問題の特徴が反映されにくい。そのため、そのような部分問題を用いて計算を行う場合、精度のよい全体解を得るため多くの反復計算を要し、計算量がより増加する可能性がある。一方、最適化装置１０では、上記のように組合せ最適化問題の特徴を反映した適切な分割が可能になるため、精度のよい全体解を得るまでの計算量の増加を抑制できる。
【0038】
（第２の実施の形態）
図３は、第２の実施の形態の最適化装置のハードウェアの一例を示す図である。
第２の実施の形態の最適化装置２０は、たとえばコンピュータであり、ＣＰＵ２１、ＲＡＭ２２、ＨＤＤ２３、画像信号処理部２４、入力信号処理部２５、媒体リーダ２６及び通信インタフェース２７を有する。上記ユニットは、バスに接続されている。
図３は、第２の実施の形態の最適化装置のハードウェアの一例を示す図である。
第２の実施の形態の最適化装置２０は、たとえばコンピュータであり、ＣＰＵ２１、ＲＡＭ２２、ＨＤＤ２３、画像信号処理部２４、入力信号処理部２５、媒体リーダ２６及び通信インタフェース２７を有する。上記ユニットは、バスに接続されている。
【0039】
ＣＰＵ２１は、プログラムの命令を実行する演算回路を含むプロセッサである。ＣＰＵ２１は、ＨＤＤ２３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ２２にロードし、プログラムを実行する。なお、ＣＰＵ２１は複数のプロセッサコアを備えてもよく、最適化装置２０は複数のプロセッサを備えてもよく、以下で説明する処理を複数のプロセッサまたはプロセッサコアを用いて並列に実行してもよい。また、複数のプロセッサの集合（マルチプロセッサ）を「プロセッサ」と呼んでもよい。
【0040】
ＲＡＭ２２は、ＣＰＵ２１が実行するプログラムやＣＰＵ２１が演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。なお、最適化装置２０は、ＲＡＭ以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。
【0041】
ＨＤＤ２３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、及び、データを記憶する不揮発性の記憶装置である。プログラムには、たとえば、最適化装置２０の制御プログラムが含まれる。なお、最適化装置２０は、フラッシュメモリやＳＳＤ（Solid State Drive）などの他の種類の記憶装置を備えてもよく、複数の不揮発性の記憶装置を備えてもよい。
【0042】
画像信号処理部２４は、ＣＰＵ２１からの命令にしたがって、最適化装置２０に接続されたディスプレイ２４ａに画像を出力する。ディスプレイ２４ａとしては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、プラズマディスプレイ（ＰＤＰ：Plasma Display Panel）、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなどを用いることができる。
【0043】
入力信号処理部２５は、最適化装置２０に接続された入力デバイス２５ａから入力信号を取得し、ＣＰＵ２１に出力する。入力デバイス２５ａとしては、マウスやタッチパネルやタッチパッドやトラックボールなどのポインティングデバイス、キーボード、リモートコントローラ、ボタンスイッチなどを用いることができる。また、最適化装置２０に、複数の種類の入力デバイスが接続されていてもよい。
【0044】
媒体リーダ２６は、記録媒体２６ａに記録されたプログラムやデータを読み取る読み取り装置である。記録媒体２６ａとして、たとえば、磁気ディスク、光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。磁気ディスクには、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤが含まれる。光ディスクには、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）が含まれる。
【0045】
媒体リーダ２６は、たとえば、記録媒体２６ａから読み取ったプログラムやデータを、ＲＡＭ２２やＨＤＤ２３などの他の記録媒体にコピーする。読み取られたプログラムは、たとえば、ＣＰＵ２１によって実行される。なお、記録媒体２６ａは、可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体２６ａやＨＤＤ２３を、コンピュータ読み取り可能な記録媒体ということがある。
【0046】
通信インタフェース２７は、ネットワーク２７ａに接続され、ネットワーク２７ａを介して他の情報処理装置と通信を行うインタフェースである。通信インタフェース２７は、スイッチなどの通信装置とケーブルで接続される有線通信インタフェースでもよいし、基地局と無線リンクで接続される無線通信インタフェースでもよい。
【0047】
図４は、第２の実施の形態の最適化装置のハードウェアの他の例を示す図である。図４において、図３に示した要素と同じ要素については同一符号が付されている。
最適化装置３０は、情報処理装置２０ａとイジングマシン２８ａを有する。情報処理装置２０ａは、インタフェース２８を有する。インタフェース２８は、イジングマシン２８ａに接続され、ＣＰＵ２１とイジングマシン２８ａとの間でデータの送受信を行う。インタフェース２８は、たとえば、ＰＣＩ（Peripheral Component Interconnect） Ｅｘｐｒｅｓｓなどの有線通信インタフェースでもよいし、無線通信インタフェースでもよい。
最適化装置３０は、情報処理装置２０ａとイジングマシン２８ａを有する。情報処理装置２０ａは、インタフェース２８を有する。インタフェース２８は、イジングマシン２８ａに接続され、ＣＰＵ２１とイジングマシン２８ａとの間でデータの送受信を行う。インタフェース２８は、たとえば、ＰＣＩ（Peripheral Component Interconnect） Ｅｘｐｒｅｓｓなどの有線通信インタフェースでもよいし、無線通信インタフェースでもよい。
【0048】
イジングマシン２８ａは、デジタル回路を用いてシミュレーテッド・アニーリング法やレプリカ交換法などを実行するハードウェアであってもよいし、量子アニーリングを行うハードウェアであってもよい。
【0049】
次に、最適化装置２０，３０の機能及び処理手順を説明する。
図５は、第２の実施の形態の最適化装置の機能例を示すブロック図である。
なお、以下では、図３に示した最適化装置２０の機能例を示すが、図４に示した最適化装置３０についても同様の機能を有する。
図５は、第２の実施の形態の最適化装置の機能例を示すブロック図である。
なお、以下では、図３に示した最適化装置２０の機能例を示すが、図４に示した最適化装置３０についても同様の機能を有する。
【0050】
最適化装置２０は、入力部３１、記憶部３２、局所解算出部３３、選択部３４、探索部３５、全体解生成部３６、出力部３７を有する。入力部３１、局所解算出部３３、選択部３４、探索部３５、全体解生成部３６、出力部３７は、たとえば、ＣＰＵ２１が実行するプログラムモジュールを用いて実装できる。なお、図４に示した最適化装置３０が用いられる場合、イジングマシン２８ａが探索部３５（またはその一部）として機能する。記憶部３２は、たとえば、ＲＡＭ２２またはＨＤＤ２３に確保した記憶領域を用いて実装できる。
【0051】
入力部３１は、たとえば、入力デバイス２５ａによって入力される選択領域情報や、たとえば、ネットワーク２７ａを介して供給される重み係数を取得する。なお、最適化装置２０が、入力された問題データに基づいて重み係数を算出するようにしてもよい。選択領域情報や重み係数は、記憶部３２に記憶される。
【0052】
局所解算出部３３は、式（１）に示したような評価関数に基づいて、記憶部３２に記憶されている重み係数を用いて、たとえば、タブーサーチ法などによって、組合せ最適化問題の局所解を算出する。なお、局所解算出部３３は、全体解生成部３６が生成した全体解を初期解として局所解を算出してもよい。
【0053】
選択部３４は、複数の状態変数のうち局所解における値が所定値（以下では１であるものとする）である状態変数と、重み係数と、入力される選択領域情報とに基づいて、複数の状態変数から部分問題の対象となる状態変数群を選択する。なお、選択部３４は、局所解の代りに、全体解生成部３６が生成した全体解に基づいて、部分問題の対象となる状態変数群を更新してもよい。
【0054】
探索部３５は、選択部３４によって選択された状態変数群についての部分問題を表す評価関数について、基底状態を探索する。探索部３５は、シミュレーテッド・アニーリング法やレプリカ交換法などのマルコフ連鎖モンテカルロ法により基底状態の探索を行ってもよいし、量子アニーリング法により基底状態の探索を行ってもよい。
【0055】
全体解生成部３６は、部分解に基づいて局所解を更新することで全体解を生成する。
出力部３７は、これまで得られた全体解のうち、評価関数の値が最小となる全体解を、最良の解として、たとえば、ディスプレイ２４ａに出力し表示させる。出力部３７は、最良の解を、記憶部３２に記憶させてもよい。
出力部３７は、これまで得られた全体解のうち、評価関数の値が最小となる全体解を、最良の解として、たとえば、ディスプレイ２４ａに出力し表示させる。出力部３７は、最良の解を、記憶部３２に記憶させてもよい。
【0056】
次に、第２の実施の形態の最適化装置２０，３０の動作（最適化装置２０，３０の制御方法）の例を説明する。
なお、以下では、計算対象の組合せ最適化問題の一例として、創薬分野で重要になるアミノ酸の安定構造探索を挙げて説明を行う。
なお、以下では、計算対象の組合せ最適化問題の一例として、創薬分野で重要になるアミノ酸の安定構造探索を挙げて説明を行う。
【0057】
図６は、組合せ最適化手法を利用したアミノ酸の安定構造の決定過程の一例を示す図である。
図６に示すような複数のアミノ酸（アミノ酸残基）が結合された直鎖構造において、以下では、各アミノ酸を１つの粒子として仮定する。図６では、８つのアミノ酸による直鎖構造における各アミノ酸が、粒子で表された例が示されている。８つの粒子において“Ｄ”はアスパラギン酸を表し、“Ａ”はアラニンを表し、“Ｙ”はチロシンを表し、“Ｑ”はグルタミンを表し、“Ｗ”はトリプトファンを表し、“Ｌ”はロイシンを表し、“Ｋ”はリシンを表す。
図６に示すような複数のアミノ酸（アミノ酸残基）が結合された直鎖構造において、以下では、各アミノ酸を１つの粒子として仮定する。図６では、８つのアミノ酸による直鎖構造における各アミノ酸が、粒子で表された例が示されている。８つの粒子において“Ｄ”はアスパラギン酸を表し、“Ａ”はアラニンを表し、“Ｙ”はチロシンを表し、“Ｑ”はグルタミンを表し、“Ｗ”はトリプトファンを表し、“Ｌ”はロイシンを表し、“Ｋ”はリシンを表す。
【0058】
アミノ酸の立体構造は、３次元格子空間の格子点にアミノ酸を配置するＬＰ（Lattice Protein）モデルを用いアミノ酸の配置位置をダイヤモンド・エンコーディング法にしたがって記述することで、イジングモデルに置き換えられる。
【0059】
図６には、イジングモデルを用いた最適化処理によって得られる安定構造の例が３次元格子空間上で示されている。最適化処理によって得られた解（安定構造）が、分子動力学計算などの詳細構造計算の際の初期解として使用される。このように、最適化処理を用いることで、アミノ酸の安定構造決定の際の総コストを低減できる。
【0060】
ダイヤモンド・エンコーディング法は、R. Babbush, and et. al., “Construction of Energy Functions for Lattice Heteropolymer Models: A Case Study in Constraint Satisfaction Programming and Adiabatic Quantum Optimization”, Advances in Chemical Physics, 155, 201-244, 13 June 2013に記載されているが、以下概略を説明する。
【0061】
図７は、ダイヤモンド・エンコーディング法の例を示す図である。
図７の例では、アミノ酸４０，４１，４２，４３，４４がそれぞれ粒子で表されている。ダイヤモンド・エンコーディング法は、直鎖構造における各アミノ酸がどの位置にあるのかをダイヤモンド格子上に埋め込んでいく手法であり、３次元の構造体が表現できる。図７では、説明を簡略化するために２次元の構造体を表現する例が示されている。各格子点は複数のアミノ酸の何れかが配置される候補である。直鎖構造の先頭のアミノ酸４０は、ダイヤモンド格子上の中心の格子点からなるシェル５０に配置され、２番目のアミノ酸４１は、シェル５０に隣接するシェル５１に属す格子点の何れかに配置される。３番目のアミノ酸４２は、シェル５１に隣接するシェル５２に属す格子点の何れかに配置され、４番目のアミノ酸４３は、シェル５２に隣接するシェル５３に属す格子点の何れかに配置される。５番目のアミノ酸４４についても、シェル５３に隣接する図示しないシェルに属す格子点の何れかに配置される。
図７の例では、アミノ酸４０，４１，４２，４３，４４がそれぞれ粒子で表されている。ダイヤモンド・エンコーディング法は、直鎖構造における各アミノ酸がどの位置にあるのかをダイヤモンド格子上に埋め込んでいく手法であり、３次元の構造体が表現できる。図７では、説明を簡略化するために２次元の構造体を表現する例が示されている。各格子点は複数のアミノ酸の何れかが配置される候補である。直鎖構造の先頭のアミノ酸４０は、ダイヤモンド格子上の中心の格子点からなるシェル５０に配置され、２番目のアミノ酸４１は、シェル５０に隣接するシェル５１に属す格子点の何れかに配置される。３番目のアミノ酸４２は、シェル５１に隣接するシェル５２に属す格子点の何れかに配置され、４番目のアミノ酸４３は、シェル５２に隣接するシェル５３に属す格子点の何れかに配置される。５番目のアミノ酸４４についても、シェル５３に隣接する図示しないシェルに属す格子点の何れかに配置される。
【0062】
アミノ酸は、１つのシェルに１つだけ配置されるという制約がある。つまり、同一シェルにおいて、２以上のアミノ酸の配置が禁止されている。なお、格子点の中には、複数のシェルにまたがって属すものがある。そのような格子点には、１つのシェルにおいてアミノ酸が配置される場合、他のシェルにおいてはアミノ酸が配置されない。
【0063】
ＬＰモデルをイジングモデルで表すため、各シェルに属す格子点には、評価関数に含まれる状態変数が割当てられる。
図８は、各シェルの格子点に割当てられた状態変数の一例を示す図である。
図８は、各シェルの格子点に割当てられた状態変数の一例を示す図である。
【0064】
シェル５０に属す格子点には、ｘ１が割当てられ、シェル５１に属す４つの格子点には、ｘ２～ｘ５が割当てられる。また、シェル５２に属す８つの格子点には、ｘ６～ｘ１３が割当てられ、シェル５３の１６個の格子点には、ｘ１４～ｘ２９が割当てられる。
【0065】
各シェルにおいて、アミノ酸が配置された格子点に割当てられた状態変数の値は、１となる。
図９は、各シェルにおけるアミノ酸の配置可能位置の一例を示す図である。
図９は、各シェルにおけるアミノ酸の配置可能位置の一例を示す図である。
【0066】
シェル５０に属す格子点にアミノ酸が配置され（ｘ１＝１）、Ｗ１２，Ｗ１３，Ｗ１４，Ｗ１５≠０である場合、シェル５１に属す各格子点にアミノ酸が配置可能である。シェル５１において、ｘ２が割当てられた格子点にアミノ酸が配置され（ｘ２＝１）、Ｗ２６，Ｗ２７，Ｗ２ １３≠０である場合、シェル５２に属し、ｘ６，ｘ７，ｘ１３が割当てられた格子点にアミノ酸が配置可能である。シェル５２において、ｘ６が割当てられた格子点にアミノ酸が配置され（ｘ６＝１）、Ｗ６ １４，Ｗ６ １８，Ｗ６ １９，Ｗ６ ２９≠０である場合、シェル５３に属し、ｘ１４，ｘ１８，ｘ１９，ｘ２９が割当てられた格子点にアミノ酸が配置可能である。
【0067】
ダイヤモンド・エンコーディング法を用いて得られるイジングモデルの評価関数は、Ｅ（ｘ）＝Ｈｏｎｅ＋Ｈｃｏｎｎ＋Ｈｏｌａｐ＋Ｈｐａｉｒ、と表せる。
Ｈｏｎｅは、１つのシェルに存在するアミノ酸が１つであるという制約条件を満たさない場合に発生する制約項であり、以下の式（２）で表せる。
Ｈｏｎｅは、１つのシェルに存在するアミノ酸が１つであるという制約条件を満たさない場合に発生する制約項であり、以下の式（２）で表せる。
【0068】
【数2】
【0069】
式（２）において、Ｎは、直鎖構造におけるアミノ酸の数（シェルの数に相当）を示し、ｘａ，ｘｂは、ｉ番目のシェルに属す格子点集合Ｑｉに含まれる２つの格子点に割当てられた状態変数である。λｏｎｅは、所定の大きさの係数である。λｏｎｅを大きく設定することで、１つのシェルにアミノ酸が２つ以上存在した場合のＨｏｎｅが大きくなる（制約が強くなる）。
【0070】
Ｈｃｏｎｎは、全てのアミノ酸はそれぞれ繋がっているという制約条件を満たさない場合に発生する制約項であり、以下の式（３）で表せる。
【0071】
【数3】
【0072】
式（３）において、ｘｄは、ｉ番目のシェルに属す格子点集合Ｑｉに含まれる格子点に割当てられた状態変数である。ｘｕは、ｉ＋１番目のシェルに属す格子点集合Ｑｉ＋１と、ｘｄが割当てられた格子点に配置されたアミノ酸に接続するアミノ酸が配置可能な格子点集合η（ｘｄ）の共通部分に含まれる格子点に割当てられた状態変数である。λｃｏｎｎは、所定の大きさの係数である。λｃｏｎｎを大きく設定することで、全てのアミノ酸はそれぞれ繋がっているという制約条件を満たさない場合のＨｃｏｎｎが大きくなる（制約が強くなる）。
【0073】
Ｈｏｌａｐは、複数のアミノ酸が、異なる格子点に配置されるという制約条件を満たさない場合に発生する制約項であり、以下の式（４）で表せる。
【0074】
【数4】
【0075】
式（４）において、νは全格子点集合Ｖに含まれる格子点を表し、ｘａ，ｘｂ∈θ（ν）は、同じ格子点νに割当てられた２つの状態変数ｘａ，ｘｂを表している。λｏｌａｐは、所定の大きさの係数である。λｏｌａｐを大きく設定することで、複数のアミノ酸が、異なる格子点に配置されるという制約条件を満たさない場合のＨｏｌａｐが大きくなる（制約が強くなる）。
【0076】
Ｈｐａｉｒは、アミノ酸同士の相互作用によって決まる最小化したいコスト項であり、以下の式（５）で表せる。
【0077】
【数5】
【0078】
式（５）において、Ｐω（ｘａ）（ｘｂ）は、ｘａ＝ｘｂ＝１のときに生じる相互作用の大きさを表す。Ｐω（ｘａ）（ｘｂ）が正の場合、ｘａとｘｂが割当てられた２つの格子点に配置されるアミノ酸間に斥力が働いていることを表し、Ｐω（ｘａ）（ｘｂ）が負の場合、ｘａとｘｂが割当てられた２つの格子点に配置されるアミノ酸間に引力が働いていることを表す。
【0079】
Ｅ（ｘ）＝Ｈｏｎｅ＋Ｈｃｏｎｎ＋Ｈｏｌａｐ＋Ｈｐａｉｒは、式（１）に示したような、重み係数を用いた表現に変形することができる。上記の各制約項などは、重み係数の値により表現される。他の重み係数よりも値が著しく大きい重み係数は、制約項が反映されている。
【0080】
ＬＰモデルに対応した部分問題の作成のための状態変数の選択指針は、以下の通りである。
単に局所解に基づいて、アミノ酸が配置されている格子点（割当てられている状態変数の値が１である格子点）だけを選択して部分問題を作成した場合、その部分問題に対して基底状態の探索を行っても、直鎖構造自体が変わらない可能性がある。そのため、現在アミノ酸が配置されていないが、配置されることで現在の直鎖構造が変化する格子点に割当てられた状態変数が優先的に選択されるようにする。これは、ある格子点に割当てられた状態変数（ｘｉ）値が１のとき、Ｗｉｊ≠０である、識別番号＝ｊの状態変数（ｘｊ）を優先的に選択することで実現できる。
単に局所解に基づいて、アミノ酸が配置されている格子点（割当てられている状態変数の値が１である格子点）だけを選択して部分問題を作成した場合、その部分問題に対して基底状態の探索を行っても、直鎖構造自体が変わらない可能性がある。そのため、現在アミノ酸が配置されていないが、配置されることで現在の直鎖構造が変化する格子点に割当てられた状態変数が優先的に選択されるようにする。これは、ある格子点に割当てられた状態変数（ｘｉ）値が１のとき、Ｗｉｊ≠０である、識別番号＝ｊの状態変数（ｘｊ）を優先的に選択することで実現できる。
【0081】
図１０は、アミノ酸が配置されている格子点の例を示す図である。
たとえば、アミノ酸が配置されている格子点６０，６１，６２に割当てられた状態変数だけが選択された場合、基底状態の探索により、格子点６０～６２の何れかに割当てられた状態変数の値が０に変わる（アミノ酸が配置されなくなる）可能性がある。しかし、立体構造自体は変わらない。
たとえば、アミノ酸が配置されている格子点６０，６１，６２に割当てられた状態変数だけが選択された場合、基底状態の探索により、格子点６０～６２の何れかに割当てられた状態変数の値が０に変わる（アミノ酸が配置されなくなる）可能性がある。しかし、立体構造自体は変わらない。
【0082】
一方、アミノ酸が配置されている格子点６３とアミノ酸が配置されていない格子点６４に割当てられた状態変数間に相互作用がある場合（重み係数が非０）、格子点６４に割当てられた状態変数を選択することで、探索により立体構造が変化する可能性がある。
【0083】
図５に示した選択部３４は、複数の状態変数のうち局所解における値が１である状態変数と重み係数に基づいて、上記のように、複数の状態変数から部分問題の対象となる状態変数群を選択する。なお、選択部３４は、選択領域情報に基づいて、状態変数群に含まれる変数の数を制限する。この点については、後述する。
【0084】
以下、第２の実施の形態の最適化装置２０，３０の制御方法の例を説明する。全体の処理の流れは、図２に示した処理の流れと同じである。以下では、第２の実施の形態の最適化装置２０，３０による状態変数群の選択処理（部分問題の作成処理）（図２のステップＳ３の処理に相当する）の例を説明する。
【0085】
図１１は、部分問題の作成処理の一例の流れを示すフローチャートである。
選択部３４は、局所解（または全体解）から、値が１の状態変数ｘｉの識別番号ｉを１つ抽出し（ステップＳ１０）、ある変数ｘをｘ＝ｉ－ｄと設定する（ステップＳ１１）。ｄは、選択領域情報に含まれる値であり、状態変数群に含まれる状態変数の数を増やすために設定される値である。ｄは、たとえば、探索部３５が処理可能な演算ビット数に応じて決定される。
選択部３４は、局所解（または全体解）から、値が１の状態変数ｘｉの識別番号ｉを１つ抽出し（ステップＳ１０）、ある変数ｘをｘ＝ｉ－ｄと設定する（ステップＳ１１）。ｄは、選択領域情報に含まれる値であり、状態変数群に含まれる状態変数の数を増やすために設定される値である。ｄは、たとえば、探索部３５が処理可能な演算ビット数に応じて決定される。
【0086】
ステップＳ１１の処理後、選択部３４は、ｘ≦ｉ＋ｄであるか否かを判定する（ステップＳ１２）。ｘ≦ｉ＋ｄであると判定した場合、選択部３４は、Ｗｘｊ≠０のうち、値が所定範囲内であるＷｘｊの識別番号ｊを、記憶部３２に保持させる（ステップＳ１３）。上記所定範囲は、たとえば、各格子点にアミノ酸が配置された際にその相互作用の影響がおよぶ格子点に割当てられた状態変数のみが選ばれるように決定される。
【0087】
たとえば、所定範囲は、外部より設定可能な任意のパラメータであるＡ～Ｄを用いて以下のような条件の何れかを満たすように決定される。
条件１：Ａ≦Ｗｘｊ≦Ｂかつ、Ｗｘｊ≠０かつＡ＝Ｂ
条件２：Ａ≦Ｗｘｊ≦Ｂかつ、Ｗｘｊ≠０かつＡ≠Ｂ
条件３：（Ａ≦Ｗｘｊ≦ＢまたはＣ≦Ｗｘｊ≦Ｄ）かつＷｘｊ≠０かつＡ≠Ｂ≠Ｃ≠Ｄ
条件４：（Ａ≦ＷｘｊまたはＷｘｊ≦Ｂ）かつＷｘｊ≠０かつＢ≦Ａ
たとえば、各アミノ酸間の相互作用によって決定される｜Ｗｘｊ｜の最大値をＪ０とすると、－Ｊ０≦Ｗｘｊ≦Ｊ０と設定される（上記の条件２でＡ＝－Ｊ０、Ｂ＝Ｊ０の場合に相当する）。なお、重み係数ごとに、異なる所定範囲が設定されていてもよい。このような所定範囲の情報は、選択領域情報に含まれる。
条件１：Ａ≦Ｗｘｊ≦Ｂかつ、Ｗｘｊ≠０かつＡ＝Ｂ
条件２：Ａ≦Ｗｘｊ≦Ｂかつ、Ｗｘｊ≠０かつＡ≠Ｂ
条件３：（Ａ≦Ｗｘｊ≦ＢまたはＣ≦Ｗｘｊ≦Ｄ）かつＷｘｊ≠０かつＡ≠Ｂ≠Ｃ≠Ｄ
条件４：（Ａ≦ＷｘｊまたはＷｘｊ≦Ｂ）かつＷｘｊ≠０かつＢ≦Ａ
たとえば、各アミノ酸間の相互作用によって決定される｜Ｗｘｊ｜の最大値をＪ０とすると、－Ｊ０≦Ｗｘｊ≦Ｊ０と設定される（上記の条件２でＡ＝－Ｊ０、Ｂ＝Ｊ０の場合に相当する）。なお、重み係数ごとに、異なる所定範囲が設定されていてもよい。このような所定範囲の情報は、選択領域情報に含まれる。
【0088】
ステップＳ１３の処理後、選択部３４は、ｘ＝ｘ＋１とし（ステップＳ１４）、ステップＳ１２からの処理を繰り返す。
図１２は、保持される識別番号ｊの一例を示す図である。
図１２は、保持される識別番号ｊの一例を示す図である。
【0089】
図１２に示されているような重み係数による係数行列において、ｘｉ＝１の場合、Ｗｉ１～Ｗｉｎのうち、値が０ではなく、－Ｊ０以上、Ｊ０以下である重み係数の２つの識別番号のうち、ｉ以外の識別番号（図１２の例ではｊ）が保持される。ｘｉ－ｄ～ｘｉ＋ｄについても値が０か１かによらず、同様の処理が行われる。
【0090】
選択部３４は、ステップＳ１２の処理で、ｘ≦ｉ＋ｄではないと判定した場合、値が１の状態変数ｘｉの識別番号ｉを全て抽出したか否かを判定する（ステップＳ１５）。全ての識別番号ｉを抽出していないと判定した場合、選択部３４は、ステップＳ１０からの処理を繰り返す。
【0091】
選択部３４は、値が１の状態変数ｘｉの識別番号ｉを全て抽出したと判定した場合、保持している識別番号ｊに基づいて、部分問題の対象となる状態変数群を決定し（ステップＳ１６）、部分問題の作成を終了する。
【0092】
たとえば、選択部３４は、保持された識別番号ｊの状態変数の全てを選択して状態変数群に含める。さらに、選択部３４は、保持されていない識別番号の状態変数を、識別番号の大きい順（または小さい順）に、探索部３５が処理可能な演算ビット数に達するまで状態変数群に追加していく。なお、選択部３４は、保持されていない識別番号の状態変数を、探索部３５が処理可能な演算ビット数に達するまで、乱数に基づいてランダムに状態変数群に追加していってもよい。
【0093】
また、選択部３４は、たとえば、探索部３５が処理可能な演算ビット数が少ない場合、保持されていない識別番号の状態変数を、部分問題の対象となる状態変数群に含めなくてもよい。その場合、選択部３４は、保持された識別番号ｊの状態変数の中から、識別番号の大きい順、小さい順、またはランダムに、探索部３５が処理可能な演算ビット数に達するまで状態変数を選択して状態変数群に含めればよい。
【0094】
また、選択部３４は、探索部３５が処理可能な演算ビット数と、保持された識別番号ｊの数との比較結果に基づいて、状態変数群に含める状態変数を決定してもよい。
図１３は、部分問題の作成処理の他の例の流れを示すフローチャートである。
図１３は、部分問題の作成処理の他の例の流れを示すフローチャートである。
【0095】
ステップＳ１５までの処理は、図１１に示した処理と同じである。ステップＳ１５の処理後、選択部３４は、保持した識別番号ｊの数が、探索部３５が処理可能な演算ビット数（ｂｉｔｍａｘ）以上であるか否かを判定する（ステップＳ１６ａ）。識別番号ｊの数が、ｂｉｔｍａｘ以上であると判定した場合、選択部３４は、ｂｉｔｍａｘ分の識別番号ｊの状態変数を、識別番号の大きい順、小さい順またはランダムに選択して状態変数群に含め（ステップＳ１６ｂ）、部分問題の作成を終了する。識別番号ｊの数が、ｂｉｔｍａｘ以上ではないと判定した場合、選択部３４は、識別番号ｊの状態変数の全ての他、保持されていない識別番号（非ｊ）の状態変数を、ｂｉｔｍａｘ分まで状態変数群に追加する。非ｊの状態変数を状態変数群に追加していく方法として、前述のように、識別番号の大きい順、小さい順、またはランダムに追加する方法がある。ステップＳ１６ｃの処理後、選択部３４は、部分問題の作成を終了する。
【0096】
（計算例１）
以下、基底状態（最適解）のときの評価関数の値（最低エネルギー）が知られている、２０個のアミノ酸（アミノ酸残基数＝２０）による直鎖構造の安定構造決定問題の計算結果を示す。問題の規模（問題のビットサイズ）は、７．８ｋビットである。最低エネルギーは、－５８２３である。
以下、基底状態（最適解）のときの評価関数の値（最低エネルギー）が知られている、２０個のアミノ酸（アミノ酸残基数＝２０）による直鎖構造の安定構造決定問題の計算結果を示す。問題の規模（問題のビットサイズ）は、７．８ｋビットである。最低エネルギーは、－５８２３である。
【0097】
使用したハードウェアは、デジタル回路を用いてシミュレーテッド・アニーリング法による基底状態の探索を行うイジングマシン２８ａ（１ｋビット規模）を含む最適化装置３０（図４）である。
【0098】
図２に示したメインループの反復回数は３回、サブループの反復回数は２回、サブループ１回当たりの、イジングマシン２８ａによる状態更新処理の反復回数は１０６、選択領域情報に含まれる前述のｄは１とした。また、部分問題の対象となる状態変数群を選択する際に用いられる重み係数の値の範囲は、たとえば、前述の条件１～条件４の何れかを満たすとともに、部分問題の対象となる状態変数群が１ｋビット以下となるようにした。
【0099】
部分問題の対象となる状態変数群を選択するときに用いる乱数のシード値をそれぞれ変えた１００回の試行について、正答率を評価した。
比較のために部分問題の対象となる全ての状態変数を乱数に基づいてランダムに選択した場合、解が局所解に落ち込んでしまい、１００回の試行に対して得られた最低の評価関数の値は、－５６２４であった。つまり正答率は０であった。なお、計算時間は、１回の試行あたり、２６．４秒であった。
比較のために部分問題の対象となる全ての状態変数を乱数に基づいてランダムに選択した場合、解が局所解に落ち込んでしまい、１００回の試行に対して得られた最低の評価関数の値は、－５６２４であった。つまり正答率は０であった。なお、計算時間は、１回の試行あたり、２６．４秒であった。
【0100】
これに対して、上記のような処理により選択部３４が部分問題を作成した場合、１００回の試行に対して８６回、評価関数の値が－５８２３となり最適解が得られた。つまり、正答率は０．８６であった。なお、計算時間は、１回の試行あたり、２６．５秒であった。
【0101】
（計算例２）
次に、最低エネルギーが知られている、４８個のアミノ酸（アミノ酸残基数＝４８）による直鎖構造の安定構造決定問題の計算結果を示す。問題の規模（問題のビットサイズ）は、３１ｋビットである。最低エネルギーは、－２６１５である。
次に、最低エネルギーが知られている、４８個のアミノ酸（アミノ酸残基数＝４８）による直鎖構造の安定構造決定問題の計算結果を示す。問題の規模（問題のビットサイズ）は、３１ｋビットである。最低エネルギーは、－２６１５である。
【0102】
使用したハードウェアは、デジタル回路を用いてシミュレーテッド・アニーリング法による基底状態の探索を行うイジングマシン２８ａ（８ｋビット規模）を含む最適化装置３０（図４）である。
【0103】
図２に示したメインループの反復回数は３回、サブループの反復回数は２回、サブループ１回当たりの、イジングマシン２８ａによる状態更新処理の反復回数は２×１０７、選択領域情報に含まれる前述のｄは１とした。また、部分問題の対象となる状態変数群を選択する際に用いられる重み係数の値の範囲は、たとえば、前述の条件１～条件４の何れかを満たすとともに、部分問題の対象となる状態変数群は８ｋビット以下となるようにした。
【0104】
部分問題の対象となる状態変数群を選択するときに用いる乱数のシード値をそれぞれ変えた１００回の試行について、正答率を評価した。
選択部３４が部分問題を作成した場合、１００回の試行に対して９９回、評価関数の値が－２６１５となり最適解が得られた。つまり、正答率は０．９９であった。なお、計算時間は、１回の試行あたり、５６．４秒であった。
選択部３４が部分問題を作成した場合、１００回の試行に対して９９回、評価関数の値が－２６１５となり最適解が得られた。つまり、正答率は０．９９であった。なお、計算時間は、１回の試行あたり、５６．４秒であった。
【0105】
図１４は、問題の分割の有無による計算量の比較例を示す図である。
図１４では、上記の２つの問題に対して、最適化装置３０が部分問題を作成して問題を解いた場合と、問題を分割せずに直接シミュレータが問題を解いた場合についての計算量の比較結果が示されている。“分割あり（Ｈａｒｄｗａｒｅ）”が、最適化装置３０が部分問題を作成して問題を解いた場合を示し、“分割なし（Ｓｉｍｕｌａｔｉｏｎ）”が、問題を分割せずに直接シミュレータが問題を解いた場合を示す。
図１４では、上記の２つの問題に対して、最適化装置３０が部分問題を作成して問題を解いた場合と、問題を分割せずに直接シミュレータが問題を解いた場合についての計算量の比較結果が示されている。“分割あり（Ｈａｒｄｗａｒｅ）”が、最適化装置３０が部分問題を作成して問題を解いた場合を示し、“分割なし（Ｓｉｍｕｌａｔｉｏｎ）”が、問題を分割せずに直接シミュレータが問題を解いた場合を示す。
【0106】
問題ビットサイズは、前述の２つの問題の規模を示している。最適化計算ビットサイズは、最適化装置３０が部分問題を作成して問題を解く場合、イジングマシン２８ａが処理可能な演算ビット数であり、問題を分割せずに直接シミュレータが問題を解く場合、問題ビットサイズと等しい。
【0107】
Ｈａｒｄｗａｒｅ呼び出し回数は、最適化装置３０が部分問題を作成して問題を解く場合、イジングマシン２８ａが使用された回数（ステップＳ４の処理回数が行われた回数）である。問題ビットサイズが３．１０Ｅ＋０４（３１ｋ）ビットの場合、Ｈａｒｄｗａｒｅ呼び出し回数は５回、問題ビットサイズが７．８０Ｅ＋０３（７．８ｋ）ビットの場合、Ｈａｒｄｗａｒｅ呼び出し回数は６回である。シミュレータが問題を解く場合のＨａｒｄｗａｒｅ呼び出し回数は１回となっている。Ｈａｒｄｗａｒｅ呼び出し１回当たりの、イジングマシン２８ａによる探索処理の反復回数は、問題ビットサイズが３．１０Ｅ＋０４ビットの場合、２．００Ｅ＋０７、問題ビットサイズが７．８０Ｅ＋０３ビットの場合、１．００Ｅ＋０６である。シミュレータが問題を解く場合の反復回数は、問題ビットサイズが３．１０Ｅ＋０４ビットの場合、７．４２Ｅ＋０６、問題ビットサイズが７．８０Ｅ＋０３ビットの場合、２．３３Ｅ＋０６である。
【0108】
正答率は、最適化装置３０が部分問題を作成して問題を解いた場合では、３．１０Ｅ＋０４ビットの問題ビットサイズの問題において０．９９、７．８０Ｅ＋０３の問題ビットサイズの問題において０．８６であった。シミュレータが問題を解いたときの正答率は、２つの規模の問題において共に０．９９であった。
【0109】
総反復回数は、最適化装置３０が部分問題を作成して問題を解いた場合では、３．１０Ｅ＋０４ビットの問題ビットサイズの問題において１．００Ｅ＋０８、７．８０Ｅ＋０３の問題ビットサイズの問題において６．００Ｅ＋０６であった。シミュレータが問題を解いたときの総反復回数は、３．１０Ｅ＋０４ビットの問題ビットサイズの問題において７．４２Ｅ＋０６、７．８０Ｅ＋０３の問題ビットサイズの問題において２．３３Ｅ＋０６であった。
【0110】
ＩＴＳ（Iteration To Solution）（９９％の確率で正解を得るまでに必要となる反復回数）は、最適化装置３０が部分問題を作成して７．８０Ｅ＋０３の問題ビットサイズの問題を解いた場合以外では、総反復回数と同じである。計算量（計算コスト）を、ＩＴＳ×最適化計算ビットサイズで表したとき、規模の大きい問題については、分割しないで解いた場合に比べて、８．００Ｅ＋１１／２．３０Ｅ＋１１、つまり３．５倍程度の計算量の増加で抑えられることがわかった。一方、規模の小さい問題については、分割しないで解いた場合に比べて、１．４１Ｅ＋１０／１．８２Ｅ＋１０、つまり０．７７倍程度となり、分割しない場合よりも計算量が減ることがわかった。
【0111】
なお、上記の説明では、計算対象の組合せ最適化問題として、複数のアミノ酸による直鎖構造の安定構造を決定する問題を例にして説明したが、巡回セールスマン問題など、他の組合せ最適化問題についても、同様に部分問題を作成して解くことができる。巡回セールスマン問題を解く場合、各状態変数は、各格子点にアミノ酸が配置されているか否かを示す代りに、ある時刻にある都市にセールスマンがいるか否かを示す。ある時刻にある都市にセールスマンがいる場合、その時刻と都市に対応した状態変数の値が、たとえば、１になる。部分問題は、図１１や図１３に示した処理と同様の処理によって作成できる。
【0112】
【0113】
選択部３４は、同一シェルに属す格子点に割当てられた状態変数の中から、状態変数群を選択してもよい。その場合、選択部３４は、値が１であるｘｉが属すシェルにおける他の状態変数の中で、Ｗｉｊ≠０で、値が所定範囲内のＷｉｊの識別番号ｊに対応するｘｊを選択する。
【0114】
ただ、識別番号ｉが、同一シェル内の状態変数の識別番号の最大値の近傍または最小値の近傍である場合、ｉ±ｄの範囲で、選択する状態変数を決めると、他のシェルの状態変数を選択してしまい制約条件違反が発生する可能性がある。
【0115】
たとえば、図８に示したようなシェル５２においてｘ６＝１の場合、ｄ＝２であるとすると、ｉ－ｄ＝６－２＝４となり、シェル５１の状態変数であるｘ４がｘ６の代りに１となる可能性がある。そこで、選択部３４は、たとえば、以下のような処理によって、状態変数を選択することで部分問題を作成する。
【0116】
図１５は、第３の実施の形態における部分問題の作成処理の一例の流れを示すフローチャートである。
まず、選択部３４は、局所解（または全体解）から、値が１の状態変数であるｘｉの識別番号ｉを１つ抽出する（ステップＳ２０）。そして、選択部３４は、保持した識別番号ｊの数が、直鎖構造に含まれるアミノ酸の数（アミノ酸残基数）よりも少ないか否かを判定する（ステップＳ２１）。
まず、選択部３４は、局所解（または全体解）から、値が１の状態変数であるｘｉの識別番号ｉを１つ抽出する（ステップＳ２０）。そして、選択部３４は、保持した識別番号ｊの数が、直鎖構造に含まれるアミノ酸の数（アミノ酸残基数）よりも少ないか否かを判定する（ステップＳ２１）。
【0117】
識別番号ｊの数が、アミノ酸残基数よりも少ないと判定した場合、選択部３４は、ｘｉが属すシェルの番号ｋを算出する（ステップＳ２２）。そして、選択部３４は、番号ｋのシェルに属す状態変数の識別番号を昇順にＫ０，Ｋ１，…，Ｋｎ－１と定義する（ステップＳ２３）。ｎは番号ｋのシェルに属す状態変数の個数である。
【0118】
たとえば、図８に示したシェル５２の例では、ｎ＝８であり、Ｋ０＝６、Ｋ１＝７、Ｋ２＝８、Ｋ３＝９、Ｋ４＝１０、Ｋ５＝１１、Ｋ６＝１２、Ｋ７＝１３である。
選択部３４は、ｉ＝Ｋｍとして、状態変数の選択範囲を決める変数をｄ（選択領域情報に含まれる）とすると、Ｋｍ－ｄ，Ｋｍ－ｄ＋１，…，Ｋｍ，…，Ｋｍ＋ｄの状態変数のうち、Ｗｉｊ≠０で、値が所定範囲内のＷｉｊの識別番号ｊを保持する（ステップＳ２４）。ここで、選択部３４は、Ｋ０～Ｋｎ－１が円順列のように繋がっているものとして、Ｋ０～Ｋｎ－１の状態変数の中からＫｍ－ｄ～Ｋｍ＋ｄの状態変数を選択候補とする。たとえば、Ｋｍ＝Ｋｎ－１の場合、Ｋｍ＋ｄは、Ｋ０＋ｄ－１となる（ｄ＝２の場合、Ｋｍ＋２＝Ｋ１となる）。Ｋｍ＝Ｋ０の場合、Ｋｍ－ｄは、Ｋｎ－１－ｄ＋１となる（ｎ＝８，ｄ＝２の場合、Ｋｍ－２＝Ｋ６となる）。
選択部３４は、ｉ＝Ｋｍとして、状態変数の選択範囲を決める変数をｄ（選択領域情報に含まれる）とすると、Ｋｍ－ｄ，Ｋｍ－ｄ＋１，…，Ｋｍ，…，Ｋｍ＋ｄの状態変数のうち、Ｗｉｊ≠０で、値が所定範囲内のＷｉｊの識別番号ｊを保持する（ステップＳ２４）。ここで、選択部３４は、Ｋ０～Ｋｎ－１が円順列のように繋がっているものとして、Ｋ０～Ｋｎ－１の状態変数の中からＫｍ－ｄ～Ｋｍ＋ｄの状態変数を選択候補とする。たとえば、Ｋｍ＝Ｋｎ－１の場合、Ｋｍ＋ｄは、Ｋ０＋ｄ－１となる（ｄ＝２の場合、Ｋｍ＋２＝Ｋ１となる）。Ｋｍ＝Ｋ０の場合、Ｋｍ－ｄは、Ｋｎ－１－ｄ＋１となる（ｎ＝８，ｄ＝２の場合、Ｋｍ－２＝Ｋ６となる）。
【0119】
図１６は、あるシェルにおいて選択候補となる状態変数の一例を示す図である。
図１６では、シェル５２において、ｘ６＝１、ｄ＝２の場合に選択候補となる状態変数の例を示している。ｘ６＝１、つまり、ｉ＝Ｋｍ＝６の場合、前述のＫｍ－ｄ～Ｋｍ＋ｄは、Ｋ６（＝１２）、Ｋ７（＝１３），Ｋ０（＝６），Ｋ１（＝７），Ｋ２（＝８）となる。つまり、ｘ１２，ｘ１３，ｘ６，ｘ７，ｘ８が選択候補となる状態変数となり、これらの識別番号（１２，１３，６，７，８）の中から、Ｗ６ｊ≠０で、値が所定範囲内のＷ６ｊの識別番号ｊが保持される。
図１６では、シェル５２において、ｘ６＝１、ｄ＝２の場合に選択候補となる状態変数の例を示している。ｘ６＝１、つまり、ｉ＝Ｋｍ＝６の場合、前述のＫｍ－ｄ～Ｋｍ＋ｄは、Ｋ６（＝１２）、Ｋ７（＝１３），Ｋ０（＝６），Ｋ１（＝７），Ｋ２（＝８）となる。つまり、ｘ１２，ｘ１３，ｘ６，ｘ７，ｘ８が選択候補となる状態変数となり、これらの識別番号（１２，１３，６，７，８）の中から、Ｗ６ｊ≠０で、値が所定範囲内のＷ６ｊの識別番号ｊが保持される。
【0120】
ステップＳ２４の処理後、ステップＳ２０からの処理が繰り返される。ステップＳ２１の処理において、保持した識別番号ｊの数がアミノ酸残基数より小さくないと判定した場合、選択部３４は、保持している識別番号ｊに基づいて、部分問題の対象となる状態変数群を決定し（ステップＳ２５）、部分問題の作成を終了する。ステップＳ２５の処理では、前述のステップＳ１６の処理と同様に、様々な決定方法を適用可能である。また、ステップＳ２５の処理の代りに、図１３に示したステップＳ１６ａ～Ｓ１６ｃの処理を適用することもできる。
【0121】
以上のような部分問題の作成手法を用いることで、ＬＰモデルとダイヤモンド・エンコーディング法を用いたアミノ酸の安定構造決定問題に適した部分問題を作成でき、さらなる計算量の抑制が期待できる。
【0122】
なお、前述のように、上記の処理内容は、最適化装置２０，３０にプログラムを実行させることで実現できる。
プログラムは、コンピュータ読み取り可能な記録媒体（たとえば、記録媒体２６ａ）に記録しておくことができる。記録媒体として、たとえば、磁気ディスク、光ディスク、光磁気ディスク、半導体メモリなどを使用できる。磁気ディスクには、ＦＤ及びＨＤＤが含まれる。光ディスクには、ＣＤ、ＣＤ－Ｒ（Recordable）／ＲＷ（Rewritable）、ＤＶＤ及びＤＶＤ－Ｒ／ＲＷが含まれる。プログラムは、可搬型の記録媒体に記録されて配布されることがある。その場合、可搬型の記録媒体から他の記録媒体（たとえば、ＨＤＤ２３）にプログラムをコピーして実行してもよい。
プログラムは、コンピュータ読み取り可能な記録媒体（たとえば、記録媒体２６ａ）に記録しておくことができる。記録媒体として、たとえば、磁気ディスク、光ディスク、光磁気ディスク、半導体メモリなどを使用できる。磁気ディスクには、ＦＤ及びＨＤＤが含まれる。光ディスクには、ＣＤ、ＣＤ－Ｒ（Recordable）／ＲＷ（Rewritable）、ＤＶＤ及びＤＶＤ－Ｒ／ＲＷが含まれる。プログラムは、可搬型の記録媒体に記録されて配布されることがある。その場合、可搬型の記録媒体から他の記録媒体（たとえば、ＨＤＤ２３）にプログラムをコピーして実行してもよい。
【0123】
以上、実施の形態に基づき、本発明の最適化装置、最適化装置の制御方法及び最適化装置の制御プログラムの一観点について説明してきたが、これらは一例にすぎず、上記の記載に限定されるものではない。
【0124】
たとえば、上記の例では、状態変数群に含まれる状態変数の数を増やすために設定される値であるｄは、固定値としたが、状態変数群に含まれる状態変数の数が、探索部３５において処理可能な演算ビット数に達するまで増加させていってもよい。
【符号の説明】
【0125】
１０ 最適化装置
１１ 記憶部
１２ 局所解算出部
１３ 選択部
１４ 探索部
１５ 全体解生成部
１１ 記憶部
１２ 局所解算出部
１３ 選択部
１４ 探索部
１５ 全体解生成部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】