以下、図面を参照しながら、本発明の実施形態について説明する。また、図面において同一の構成要素は、同じ番号を付し、説明は、適宜省略する。
(システム構成)
図1は、情報処理システム100の構成例を示したブロック図である。図1の情報処理システム100は、管理サーバ1と、ネットワーク2と、複数の計算サーバ(情報処理装置)3(3a~3c)と、複数のケーブル4(4a~4c)と、スイッチ5を備えている。また、図1には、情報処理システム100と通信可能な情報端末6が示されている。管理サーバ1、複数の計算サーバ3(3a~3c)、情報端末6は、ネットワーク2を介して互いにデータ通信をすることができる。ネットワーク2は、例えば、複数のコンピュータネットワークが相互に接続されたインターネットである。ネットワーク2は、通信媒体として有線、無線、または、これらの組み合わせを用いることができる。また、ネットワーク2で使われる通信プロトコルの例としては、TCP/IPがあるが、通信プロトコルの種類については特に問わない。
また、複数の計算サーバ3(3a~3c)は、それぞれケーブル4(4a~4c)を介してスイッチ5に接続されている。複数のケーブル4(4a~4c)およびスイッチ5は、計算サーバ間のインターコネクトを形成している。複数の計算サーバ3(3a~3c)のそれぞれは、当該インターコネクトを介して互いにデータ通信をすることも可能である。スイッチ5は、例えば、Infinibandのスイッチであり、ケーブル4a~4cは、例えば、Infinibandのケーブルである。ただし、Infinibandのスイッチ/ケーブルの代わりに、有線LANのスイッチ/ケーブルを使ってもよい。ケーブル4a~4cおよびスイッチ5で使われる通信規格および通信プロトコルについては、特に問わない。情報端末6の例としては、ノートPC、デスクトップPC、スマートフォン、タブレット、車載端末などが挙げられる。
組合せ最適化問題の求解では、並列的な処理および/または処理の分散化を行うことができる。したがって、複数の計算サーバ3(3a~3c)のそれぞれおよび/または計算サーバ3(3a~3c)のプロセッサは、一部の計算処理の一部のステップを分担して実行してもよいし、異なる変数について同様の計算処理を並列的に実行してもよい。管理サーバ1は、例えば、ユーザによって入力された組合せ最適化問題を各計算サーバ3に処理可能な形式に変換し、計算サーバ3を制御する。そして、管理サーバ1は、各計算サーバ3から計算結果を取得し、集約した計算結果を組合せ最適化問題の解に変換する。これにより、ユーザは、組合せ最適化問題の解を得ることができる。組合せ最適化問題の解は、最適解と、最適解に近い近似解とを含むものとする。
図1には、3台の計算サーバ3(3a~3c)が示されているが、これは、情報処理システム100に含まれる計算サーバ3の台数を限定することを意図していない。また、組合せ最適化問題の求解に使われる計算サーバ3の台数についても特に問わない。例えば、情報処理システム100に含まれる計算サーバ3は1台であってもよい。また、情報処理システム100に含まれる複数の計算サーバ3のうち、1台の計算サーバ3を使って組合せ最適化問題の求解を行ってもよい。また、情報処理システム100に、数百台以上の計算サーバ3が含まれていてもよい。計算サーバ3は、データセンターに設置されたサーバであってもよいし、オフィスに設置されたデスクトップPCであってもよい。また、計算サーバ3は異なるローケーションに設置された複数の種類のコンピュータであってもよい。計算サーバ3として使われる情報処理装置の種類については特に問わない。例えば、計算サーバ3は、汎用的なコンピュータであってもよいし、専用の電子回路または、これらの組合せであってもよい。
図2は、管理サーバ1の構成例を示したブロック図である。図2の管理サーバ1は、例えば、中央演算処理装置(CPU)とメモリとを含むコンピュータである。管理サーバ1は、プロセッサ10と、記憶部14と、通信回路15と、入力回路16と、出力回路17とを備えている。プロセッサ10、記憶部14、通信回路15、入力回路16、出力回路17は、互いにバス20を介して接続されているものとする。プロセッサ10は、内部の機能構成として、管理部11と、変換部12と、制御部13を含んでいる。
プロセッサ10は、演算を実行し、管理サーバ1の制御を行う電子回路である。プロセッサ10として、例えば、CPU、マイクロプロセッサ、ASIC、FPGA、PLDまたはこれらの組み合わせを用いることができる。管理部11は、ユーザの情報端末6を介して管理サーバ1の操作を行うためのインタフェースを提供する。管理部11が提供するインタフェースの例としては、API、CLI、ウェブページなどが挙げられる。例えば、ユーザは、管理部11を介して組合せ最適化問題の情報の入力を行ったり、計算された組合せ最適化問題の解の閲覧および/またはダウンロードを行ったりすることができる。変換部12は、組合せ最適化問題に関するパラメータを入力して、入力したパラメータを各計算サーバ3が処理可能な形式に変換する。制御部13は、各計算サーバ3に制御指令を送信する。制御部13が各計算サーバ3から計算結果を取得した後、変換部12は、複数の計算結果を集約し、組合せ最適化問題の解に変換し、組合せ最適化問題の解を出力する。
記憶部14は、管理サーバ1のプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、OSとアプリケーションの両方を含むものとする。記憶部14は、揮発性メモリ、不揮発性メモリ、またはこれらの組み合わせであってもよい。揮発性メモリの例としては、DRAM、SRAMなどがある。不揮発性メモリの例としては、NANDフラッシュメモリ、NORフラッシュメモリ、ReRAM、MRAMが挙げられる。また、記憶部14として、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置を使ってもよい。
通信回路15は、ネットワーク2に接続された各装置との間でデータの送受信を行う。通信回路15は、例えば、有線LANのNIC(Network Interface Card)である。ただし、通信回路15は、無線LANなど、その他の種類の通信回路であってもよい。入力回路16は、管理サーバ1へのデータ入力を実現する。入力回路16は、外部ポートとして、例えば、USB、PCI-Expressなどを備えているものとする。図2の例では、操作装置18が入力回路16に接続されている。操作装置18は、管理サーバ1に情報を入力するための装置である。操作装置18は、例えば、キーボード、マウス、タッチパネル、音声認識装置などであるが、これに限られない。出力回路17は、管理サーバ1からのデータ出力を実現する。出力回路17は、外部ポートとしてHDMI(登録商標)、DisplayPortなどを備えているものとする。図2の例では、表示装置19が出力回路17に接続されている。表示装置19の例としては、LCD(液晶ディスプレイ)、有機EL(有機エレクトロルミネッセンス)ディスプレイ、プロジェクタがあるが、これに限られない。
管理サーバ1の管理者は、操作装置18および表示装置19を使って、管理サーバ1のメンテナンスを行うことができる。なお、操作装置18および表示装置19は、管理サーバ1に組み込まれたものであってもよい。また、必ず管理サーバ1に操作装置18および表示装置19が接続されていなくてもよい。例えば、管理者は、ネットワーク2と通信可能な情報端末を用いて管理サーバ1のメンテナンスを行ってもよい。
図3は、管理サーバ1の記憶部14に保存されるデータの例を示している。図3の記憶部14には、問題データ14Aと、計算データ14Bと、管理プログラム14Cと、変換プログラム14Dと、制御プログラム14Eとが保存されている。例えば、問題データ14Aは、組合せ最適化問題のデータを含む。例えば、計算データ14Bは、各計算サーバ3から収集された計算結果を含む。例えば、管理プログラム14Cは、上述の管理部11の機能を実現するプログラムである。例えば、変換プログラム14Dは、上述の変換部12の機能を実現するプログラムである。例えば、制御プログラム14Eは、上述の制御部13の機能を実現するプログラムである。
図4は、計算サーバ3aの構成例を示したブロックである。図4には、例示的に計算サーバ3aの構成が示されている。他の計算サーバ3は、計算サーバ3aと同様の構成であってもよいし、計算サーバ3aと異なる構成であってもよい。計算サーバ3aは、例えば、第1ベクトルと、第2ベクトルと、第3ベクトルの計算を単独で、または、他の計算サーバ3と分担して実行する情報処理装置である。計算サーバ3aは、さらに第1ベクトルのそれぞれの要素を符号関数で変換した第4ベクトルの計算を行ってもよい。第3ベクトルのそれぞれの要素の値は、例えば、イジングモデルのエネルギー式から導出された式によって求められる。例えば、イジングモデルのエネルギー式を、すべての項に含まれる変数について偏微分した形式の式(基本式とよぶ)に基づいて第3ベクトルのそれぞれの要素を計算することができる。
ここで、第1ベクトルは、変数xi(i=1、2、・・・、N)を要素とするベクトルである。また、第2ベクトルは、変数yi(i=1、2、・・・、N)を要素とするベクトルである。第3ベクトルは、変数zi(i=1、2、・・・、N)を要素とするベクトルである。第4ベクトルは、第1ベクトルの要素を第1値もしくは第1値より大きい第2値のいずれかの値をとる第1関数で変換したベクトルである。上述の符号関数は、第1関数の一例である。なお、変数xi、yi、ziの詳細については、後述する。
計算サーバ3aは、例えば、通信回路31と、共有メモリ32と、プロセッサ33A~33Dと、ストレージ34と、ホストバスアダプタ35とを備えている。通信回路31、共有メモリ32、プロセッサ33A~33D、ストレージ34、ホストバスアダプタ35は、バス36を介して互いに接続されているものとする。
通信回路31は、ネットワーク2に接続された各装置との間でデータの送受信を行う。通信回路31は、例えば、有線LANのNIC(Network Interface Card)である。ただし、通信回路31は、無線LANなど、その他の種類の通信回路であってもよい。共有メモリ32は、プロセッサ33A~33Dからアクセス可能なメモリである。共有メモリ32の例としては、DRAM、SRAMなどの揮発性メモリが挙げられる。ただし、共有メモリ32として、不揮発性メモリなどその他の種類のメモリが使われてもよい。プロセッサ33A~33Dは、共有メモリ32を介してデータの共有を行うことができる。なお、必ず計算サーバ3aのすべてのメモリが共有メモリとして構成されていなくてもよい。例えば、計算サーバ3aの一部のメモリは、いずれかのプロセッサのみからアクセスできるローカルメモリとして構成されていてもよい。
プロセッサ33A~33Dは、計算処理を実行する電子回路である。プロセッサは、例えば、CPU(Central Processing Unit)、GPU(Graphics Processing Unit)、FPGA(Field-Programmable Gate Array)、ASIC(Application Specific Integrated Circuit)のいずれであってもよいし、これらの組合せであってもよい。また、プロセッサは、CPUコアまたはCPUスレッドであってもよい。プロセッサがCPUである場合、計算サーバ3aが備えるソケット数については、特に問わない。また、プロセッサは、PCI expressなどのバスを介して計算サーバ3aのその他の構成要素に接続されていてもよい。
図4の例では、計算サーバ3aが4つのプロセッサを備えている。ただし、1台の計算サーバ3aが備えているプロセッサの数はこれとは異なっていてもよい。例えば、計算サーバ3aに実装されているプロセッサの数および/または種類が異なっていてもよい。
作用演算部51は、解きたい組合せ最適化問題の目的関数を、すべての項に含まれる変数について偏微分した形式の基本式に基づき、第3ベクトルのそれぞれの要素を更新するように構成されている。ここで、基本式の変数は、第1ベクトルの要素、または、第1ベクトルの要素を第1値もしくは第1値より大きい第2値のいずれかの値をとる第1関数で変換した第4ベクトルの要素である。更新部50は、例えば、第1ベクトルの要素に、第2ベクトルの対応する要素、または、第2ベクトルの対応する要素に重み付けした値を加算することによって第1ベクトルの要素を更新し、値が第1値より小さい第1ベクトルの要素を第1値以上しきい値以下の何れかの値に変更し、値が第2値より大きい第1ベクトルの要素をしきい値以上第2値以下の値に変更し、第2ベクトルの要素に、更新回数に応じて単調増加または単調減少する第1係数と第1ベクトルの対応する要素との積を重み付けした値、ならびに、第3ベクトルの対応する要素を重み付けした値を加算することによって第2ベクトルの要素を更新するように構成されている。なお、しきい値は、第1値と第2値との間の値である。目的関数として、例えば、イジングモデルのエネルギー式を使うことができる。ここで、イジングモデルは、多体相互作用を有するものであってもよい。また、第1値として-1、第2値として+1、しきい値として0を使うことができる。ただし、しきい値、第1値および/または第2値はその他の値であってもよい。
図4の例では、プロセッサ33A~33Cが更新部50に相当しており、プロセッサ33Dが作用演算部51に相当している。ただし、図4に示した更新部50/作用演算部51とプロセッサとの対応関係は、一例にしかすぎない。したがって、更新部50/作用演算部51とプロセッサとの対応関係はこれとは異なっていてもよい。また、更新部50/作用演算部51に割り当てられるプロセッサ数については、特に限定しない。後述するように、同一のプロセッサが更新部50および作用演算部51の役割を兼ね備えていてもよい。計算サーバ3aに複数の種類のプロセッサ(例えば、CPU、GPU、FPGA)が実装されている場合には、異なる種類のプロセッサを更新部50および作用演算部51に割り当ててもよい。
ストレージ34は、計算サーバ3aのプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、OSとアプリケーションの両方を含むものとする。ストレージ34は、揮発性メモリ、不揮発性メモリ、またはこれらの組み合わせであってもよい。揮発性メモリの例としては、DRAM、SRAMなどがある。不揮発性メモリの例としては、NANDフラッシュメモリ、NORフラッシュメモリ、ReRAM、MRAMが挙げられる。また、ストレージ34として、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置が使われてもよい。
ホストバスアダプタ35は、計算サーバ3間のデータ通信を実現する。ホストバスアダプタ35は、ケーブル4aを介してスイッチ5に接続されている。ホストバスアダプタ35は、例えば、HCA(Host Channel Adaptor)である。ホストバスアダプタ35、ケーブル4a、スイッチ5で高スループットを実現可能なインターコネクトを形成することにより、並列的な計算処理の速度を向上させることができる。
図5は、計算サーバ3のストレージに保存されるデータの例を示している。図5のストレージ34には、計算データ34Aと、計算プログラム34Bと、制御プログラム34Cとが保存されている。計算データ34Aは、計算サーバ3aの計算途中のデータまたは計算結果を含んでいる。なお、計算データ34Aの少なくとも一部は、共有メモリ32、プロセッサのキャッシュ、プロセッサのレジスタなど、異なる記憶階層に保存されていてもよい。計算プログラム34Bは、所定のアルゴリズムに基づき、各プロセッサにおける計算処理および、共有メモリ32およびストレージ34へのデータの保存処理を実現するプログラムである。制御プログラム34Cは、管理サーバ1の制御部13から送信された指令に基づき、計算サーバ3aを制御し、計算サーバ3aの計算結果を管理サーバ1に送信するプログラムである。
(組合せ最適化問題)
次に、組合せ最適化問題の求解に関連する技術について説明する。組合せ最適化問題を解くために使われる情報処理装置の一例として、イジングマシンが挙げられる。イジングマシンとは、イジングモデルの基底状態のエネルギーを計算する情報処理装置のことをいう。これまで、イジングモデルは、主に強磁性体や相転移現象のモデルとして使われることが多かった。しかし、近年、イジングモデルは、組み合わせ最適化問題を解くためのモデルとしての利用が増えている。下記の式(1)は、イジングモデルのエネルギーを示している。
ここで、s
i、s
jはスピンである、スピンは、+1または-1のいずれかの値をとる2値変数である。Nは、スピンの数である。h
iは、各スピンに作用する局所磁場である。Jは、スピン間における結合係数の行列である。行列Jは、対角成分が0である実対称行列となっている。したがって、J
ijは行列Jのi行j列の要素を示している。なお、式(1)のイジングモデルは、スピンについての2次式となっているが、スピンの3次以上の項を含む拡張されたイジングモデル(多体相互作用を有するイジングモデル)を使ってもよい。多体相互作用を有するイジングモデルについては、後述する。
式(1)のイジングモデルを使うと、エネルギーEIsingを目的関数とし、エネルギーEIsingを可能な限り小さくする解を計算することができる。イジングモデルの解は、スピンのベクトル(s1、s2、・・・、sN)の形式で表される。特に、エネルギーEIsingが最小値となるベクトル(s1、s2、・・・、sN)は、最適解とよばれる。ただし、計算されるイジングモデルの解は、必ず厳密な最適解でなくてもよい。以降では、イジングモデルを使ってエネルギーEIsingが可能な限り小さくなる近似解(すなわち、目的関数の値が可能な限り最適値に近くなる近似解)を求める問題をイジング問題とよぶものとする。
式(1)のsiは、スピンを表す2値変数であるため、式(1+si)/2を使うことにより、組合せ最適化問題で使われる離散変数(ビット)との変換を容易に行うことができる。したがって、組合せ最適化問題をイジング問題に変換し、イジングマシンに計算を行わせることにより、組合せ最適化問題の解を求めることが可能である。0または1のいずれかの値をとる離散変数(ビット)を変数とする2次の目的関数を最小化する解を求める問題は、QUBO(Quadratic Unconstrained Binary Optimization、制約なし2値変数2次最適化)問題とよばれる。式(1)で表されるイジング問題は、QUBO問題と等価であるといえる。
例えば、量子アニーラ、コヒーレントイジングマシン、量子分岐マシンなどがイジングマシンのハードウェア実装として提案されている。量子アニーラは、超伝導回路を使って量子アニーリングを実現する。コヒーレントイジングマシンは、光パラメトリック発振器で形成されたネットワークの発振現象を利用する。量子分岐マシンは、カー効果を有するパラメトリック発振器のネットワークにおける量子力学的な分岐現象を利用する。これらのハードウェア実装は、計算時間の大幅な短縮を実現する可能性がある一方、大規模化や安定的な運用が難しいという課題もある。
そこで、広く普及しているデジタルコンピュータを使ってイジング問題の求解を行うことも可能である。デジタルコンピュータは、上述の物理的現象を使ったハードウェア実装と比べ、大規模化と安定運用が容易である。デジタルコンピュータでイジング問題の求解を行うためのアルゴリズムの一例として、シミュレーティッドアニーリング(SA)が挙げられる。シミュレーティッドアニーリングをより高速に実行する技術の開発が行われている。ただし、一般のシミュレーティッドアニーリングはそれぞれの変数が逐次更新される逐次更新アルゴリズムであるため、並列化による計算処理の高速化は難しい。
(シミュレーテッド分岐アルゴリズム)
上述の課題を踏まえ、デジタルコンピュータにおける並列的な計算によって、規模の大きい組合せ最適化問題の求解を高速に行うことが可能なシミュレーテッド分岐アルゴリズムが提案されている(例えば、上述の非特許文献7)。以降では、シミュレーテッド分岐アルゴリズムを使って組合せ最適化問題を解く情報処理装置および電子回路について説明する。
はじめに、シミュレーテッド分岐アルゴリズムの概要について述べる。シミュレーテッド分岐アルゴリズムでは、それぞれがN個の要素に対応する2つの変数x
iおよび変数y
iを用いる。なお、変数x
iを第1変数、変数y
iを第2変数と呼ぶ場合もある。ここで、シミュレーテッド分岐アルゴリズムにおいて、N個の要素のそれぞれは、仮想的な粒子を表している。N個の要素は、最適化問題を表すイジングモデルにおけるN個のスピンに対応する。第1変数x
iは、N個の粒子のうちのi番目の粒子の位置を表す。第2変数y
iは、i番目の粒子の運動量を表す。Nは、イジングモデルに含まれるスピンの数を表し、2以上の整数である。iは、1以上、N以下の任意の整数を表し、スピンを特定するインデックスを表す。そして、シミュレーテッド分岐アルゴリズムでは、それぞれN個ある2つの変数x
i,y
i(i=1、2、・・・、N)について、下記の式(2)の連立常微分方程式を数値的に解く。変数x
i,y
iは、いずれも、実数で表される連続変数である。
ここで、Hは、下記の式(3)のハミルトニアンである。係数Dは、予め定められた定数であり、離調(detuning)に相当する。係数p(t)は、ポンピング振幅(pumping amplitude)に相当し、シミュレーテッド分岐アルゴリズムの計算時に更新回数に応じて値が単調増加する。tは、時刻を表す変数である。係数p(t)の初期値は0に設定されていてもよい。係数p(t)は、第1係数に相当する。係数Kは、正のカー係数(Kerr coefficient)に相当する。f
iは、下記の式(4)で表される外力である。式(4)のz
iは、式(3)の中のエネルギーE
Isingに対応する項の括弧内を変数x
iで偏微分した式となっている。
ここで、係数cとして、定数係数を使うことができる。この場合、係数cの値を、シミュレーテッド分岐アルゴリズムによる計算を実行する前に決める必要がある。例えば、計算の精度を得るために、係数cをJ
(2)行列の最大固有値の逆数に近い値に設定することができる。例えば、c=0.5D√(N/2n)という値を使うことができる。ここで、nは、組合せ最適化問題に係るグラフのエッジ数である。また、α(t)は、p(t)とともに増加する係数である。例えば、α(t)として、√(p(t))を使うことができる。
なお、シミュレーテッド分岐アルゴリズムを使うことにより、3次以上の目的関数を有する組合せ最適化問題を解くことも可能である。2値変数を変数とする3次以上の目的関数を最小化する変数の組合せを求める問題は、HOBO(Higher Order Binary Optimization)問題とよばれる。HOBO問題を扱う場合、高次へ拡張されたイジングモデルにおけるエネルギー式として、下記の式(5)を使うことができる。
ここで、J
(n)はn階テンソルであり、式(1)の局所磁場h
iと結合係数の行列Jを一般化させたものである。例えば、テンソルJ
(1)は、局所磁場h
iのベクトル(第6ベクトルとよぶ)に相当する。n階テンソルJ
(n)では、複数の添え字に同じ値があるとき、要素の値は0となる。式(5)では、3次の項までが示されているが、それより高次の項も式(5)と同様に定義することができる。式(5)は多体相互作用を含むイジングモデルのエネルギーに相当している。
なお、QUBOと、HOBOはいずれも、制約なし多項式2値変数最適化(PUBO:Polynomial Unconstrained Binary Optimization)の1種であるといえる。すなわち、PUBOのうち、2次の目的関数を有する組合せ最適化問題は、QUBOである。また、PUBOのうち、3次以上の目的関数を有する組合せ最適化問題は、HOBOであるといえる。
シミュレーテッド分岐アルゴリズムを使ってHOBO問題を解く場合、上述の式(3)のハミルトニアンHを下記の式(6)に、上述の式(4)の外力f
iを下記の式(7)にそれぞれ置き換えればよい。
例えば、(7)の2番目の式z
iを用いて、第3ベクトルのそれぞれの要素を計算することができる。この式は、(6)の2番目の式を、すべての項に含まれる変数x
iについて偏微分した形式をとっている。また、第1ベクトルの要素を変数としている。このように、ハミルトニアンが多体相互作用(3階以上のテンソル)の項を含んでもよい。また、ハミルトニアンとして、多体相互作用(3階以上のテンソル)の項を含まないものを使ってもよい。(7)の2番目の式z
iは、ハミルトニアンの中のイジングモデルのエネルギーに対応する項から導出された基本式の一例である。すなわち、第1値は-1、第2値は1であってもよく、目的関数は、イジングモデルのエネルギー式に相当する項を含んでいてもよい。この場合、目的関数は、多体相互作用の項を含んでいてもよい。
シミュレーテッド分岐アルゴリズムでは、p(t)の値を初期値(例えば、0)から所定の値まで増加させた後における変数xiの符号に基づき、スピンsiの値を求めることができる。例えば、xi>0のときsgn(xi)=1、xi<0のときsgn(xi)=-1となる符号関数を使うと、p(t)の値が所定の値まで増加したとき、変数xiを符号関数で変換することによってスピンsiの値を求めることができる。符号関数として、例えば、xi≠0のときに、sgn(xi)=xi/|xi|、xi=0のときにsgn(xi)=1または-1になる関数を使うことができる。すなわち、更新部50は、値が第1値と第2値の間にある第3値より小さい第1ベクトルの要素を第1値に変換し、値が第3値より大きい第1ベクトルの要素を第2値に変換することによって、組合せ最適化問題の解を求めるように構成されていてもよい。例えば、更新部50は、正値である第1ベクトルの要素を+1に変換し、負値である第1ベクトルを-1に変換することによって、組合せ最適化問題の解を求めるように構成されていてもよい。更新部50が組合せ最適化問題の解(例えば、イジングモデルのスピンsi)を求めるタイミングについては、特に問わない。例えば、更新部50は、第1ベクトル、第2ベクトル、第3ベクトルの更新回数または第1係数pの値がしきい値より大きくなったときに組合せ最適化問題の解を求めるように構成されていてもよい。イジング問題を解く場合、組合せ最適化問題の解は、イジングモデルのスピンsiに相当する。
(シミュレーテッド分岐アルゴリズムの演算)
例えば、シンプレクティック・オイラー法を使うと、式(2)、(3)、(4)または、式(2)、(6)、(7)によって与えられる微分方程式を解くことができる。下記の式(8)に示されているように、シンプレクティック・オイラー法を使う場合、微分方程式が離散的な漸化式に書き換えられる。
ここで、tは、時刻であり、Δtは、時間ステップ(単位時間、時間刻み幅)である。式(8)の非線形項Kx
2
i(t+Δt)は、計算中に変数x
iが発散するのを防止する。
計算サーバ3では、式(8)のアルゴリズムに基づき、それぞれN個ある2つの変数xi,yi(i=1、2、・・・、N)を更新してもよい。すなわち、計算サーバ3が更新するデータには、変数xi(i=1、2、・・・、N)を要素とする第1ベクトル(x1,x2,・・・,xN)と、変数yi(i=1、2、・・・、N)を要素とする第2ベクトル(y1,y2,・・・,yN)と、変数zi(i=1、2、・・・、N)を要素とする第3ベクトル(z1,z2,・・・,zN)が含まれていてもよい。計算サーバ3は、式(8)のアルゴリズムに基づき、第3ベクトルのそれぞれの要素zi(i=1、2、・・・、N)と、第1のベクトルのそれぞれの要素xi(i=1、2、・・・、N)と、第2ベクトルのそれぞれの要素yi(i=1、2、・・・、N)を更新することができる。
式(8)を参照すると、外力項fiに含まれる、行列またはテンソルの積和演算以外は、1種類の添え字(i)しか現れていないことがわかる。このため、式(8)のうち、1種類の添え字(i)しか現れていない部分の演算を並列化することによって、計算時間を短縮することができる。
なお、式(8)では、微分方程式との対応関係を示すために、時刻tおよび時間ステップΔtが使われている。ただし、実際にシンプレクティック・オイラー法をソフトウェアまたはハードウェアに実装する際は、必ず明示的なパラメータとして時刻tおよび時間ステップΔtが含まれていなくてもよい。例えば、時間ステップΔtを1とすれば、実装時のアルゴリズムから時間ステップΔtを除去することが可能である。アルゴリズムを実装する際に、明示的なパラメータとして時刻tを含めない場合には、式(8)において、xi(t+Δt)をxi(t)の更新後の値であると解釈すればよい。すなわち、上述の式(8)および以降の各式における“t”は、更新前の変数の値、“t+Δt”は、更新後の変数の値を示すものとする。
次に、シミュレーテッド分岐アルゴリズムをシンプレクティック・オイラー法によってデジタルコンピュータに実装し、組合せ最適化問題を解いたときの結果を示す。以降では、最大カット問題のベンチマークセット(G-set)のG22を1000回解いた場合におけるカット数の平均値と最大値を示す。最大カット問題とは、分割時にカットされるエッジの重みの合計値を最大化するよう、重み付きグラフのノードを2つのグループに分割する問題である。最大カット問題は、組合せ最適化問題の一種である。
図6Aおよび図6Bは、上述の式(8)のアルゴリズムを使ったときにおける結果を示している。時間ステップをΔt=0.5、合計時間ステップ数を100、1000、10000、100000として計算を行っている。係数について、D=K=1、c=0.5D√(N/2n)を使った。なお、nにはG22のグラフのエッジ数、19990が代入される。時間ステップ数の増加に応じて係数p(第1係数)の値を0から1に線形に増加させた。また、変数xiの初期値に0を設定し、変数yiの初期値として[-0.1,0.1]の範囲の擬似乱数を設定した。
図6Aは、カット数の平均値を示している。一方、図6Bは、カット数の最大値を示している。図6Aおよび図6Bのいずれのグラフにおいても、縦軸はカット数に、横軸は時間ステップ数にそれぞれ対応している。図6Aおよび図6Bの両グラフにある水平方向の破線Cmaxは、G22で知られている最大カット数13359を示している。カット数が破線Cmaxに近いほど、最適解に近い結果が得られているといえる。図6Aおよび図6Bを参照すると、合計時間ステップ数が大きくなっても最大カット数に到達していない。
(アルゴリズムの改良)
図7は、式(8)のアルゴリズムの分岐現象を表す図である。シミュレーテッド分岐アルゴリズムをシンプレクティック・オイラー法によって解いた場合、系のパラメータの変化に伴い、安定運動状態が1個から複数個へと分岐する。式(8)のアルゴリズムの分岐現象では、変数xiが1より大きい領域、または、-1より小さい領域まで広がっている。
図8は、改良したアルゴリズムの分岐現象を表す図である。そこで、以下のように、式(8)のアルゴリズムを改良した。具体的には、図8に示すように、更新によって変数xiの絶対値が1より大きくなった場合、変数xiは、符号が変更されずに、絶対値を0以上1以下の値に変更される。例えば、更新によってxi>1となった場合、変数xiの値はwに設定される。ここで、wは、0以上1以下の値である。また、更新によってxi<-1となった場合、変数xiの値は-wに設定される。これにより、変数xiの符号を保ちつつ、変数xiを常に-1以上、1以下の範囲に保つことができる。
例えば、wは、0以上1以下の予め定められた値であってもよい。また、wは、0以上1以下の範囲内の所定区間[w1,w2]で一様な確率で発生する乱数に応じた値であってもよい。なお、w1は、0以上w2以下となり、w2は、w1以上1以下となる。
また、wは、N個の変数xiの大きさの平均を表す指標値であってもよい。N個の変数xiの大きさの平均を表す指標値は、例えば、直前のN個の変数xiの二乗平均平方根または絶対値平均である。また、例えば、wは、N個の変数xiの大きさの平均を表す指標値以上且つ1以下の、乱数により定まる値であってもよい。
また、例えば、wは、0から1以下の間で初期時刻から終了時刻まで時間経過に従って増加する増加係数であってもよい。増加係数は、例えば、初期時刻で0、終了時刻で1となる時刻を変数とする一次関数、または、一次関数の平方根である。また、例えば、wは、増加係数以上、1以下の、乱数により定まる値であってもよい。
さらに、更新によってxi>1となった場合、変数xiに対応する変数yiに係数rを乗算してもよい。すなわち、更新部50は、値が第1値より小さい第1ベクトルの要素に対応する第2ベクトルの要素、または、第2値より大きい第1ベクトルの要素に対応する第2ベクトルの要素を、もとの前記第2ベクトルの要素に、第2係数を乗じた値に更新するように構成されていてもよい。例えば、更新部50は、値が-1より小さい第1ベクトルの要素に対応する第2ベクトルの要素、または、値が1より大きい第1ベクトルの要素に対応する第2ベクトルの要素を、もとの第2ベクトルの要素に第2係数を乗じた値に更新するように構成されていてもよい。ここで、第2係数は上述の係数rに相当する。
また、更新によって変数xiの絶対値が1より大きくなった場合、変数yiは、0または予め定められた値に変更されてもよい。また、更新によってxi>1となった場合、変数xiに対応する変数yiの値を擬似乱数に設定してもよい。例えば、[-0.1,0.1]の範囲の乱数を使うことができる。すなわち、更新部50は、値が第1値より小さい第1ベクトルの要素に対応する第2ベクトルの要素の値、または、値が第2値より大きい第1ベクトルの要素に対応する第2ベクトルの要素の値を、擬似乱数に設定するように構成されていてもよい。
(第1アルゴリズム)
以上のようにして|x
i|>1とならないように更新すると、式(8)の非線形項Kx
2
i(t+Δt)を除去しても、x
iの値が発散することはなくなる。したがって、式(8)のアルゴリズムに代わり、下記の式(9)の第1アルゴリズムを使うことが可能となる。
なお、QUBO問題の場合、式(9)のz
i(t+Δt)は、下記の式(10)により表される。
上述の式(9)の第1アルゴリズムでは、必ずしも擬似乱数を使う必要がない。また、式(9)の第1アルゴリズムは、式(8)と同様、ハミルトン方程式を解くものであり、変数yiは運動量に相当する。そのため、シンプレクティック・オイラー法を使うことにより、時間ステップΔtとして小さな値を使わなくても、安定的に解を求めることができる。また、式(9)の第1アルゴリズムにおいても、3次以上の目的関数を有する組合せ最適化問題を解くことが可能である。
図9Aおよび図9Bは、式(9)の第1アルゴリズムを使ってG-setのG22を1000回解いたときの結果を示している。図9Aおよび図9Bでは、式(9)の第1アルゴリズムを使っている。第2係数rの値は0に設定され、時間ステップはΔt=1に設定されている。wは、区間[xav2,1]における一様乱数とした。xav2は、N個の変数xiの二乗平均平方根である。その他の計算条件は、図6Aおよび図6Bと同様であるものとする。なお、式(9)では、非線形項がないため、時間ステップΔtを図6Aおよび図6Bの倍に設定することができた。
図9Aは、カット数の平均値を示している。一方、図9Bは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Cmaxの定義は、図6Aおよび図6Bと同様である。また、図9Aおよび図9Bの両グラフにおいて実線で示されたデータは、式(9)の第1アルゴリズムの適用時の結果に相当している。一方、図9Aおよび図9Bの両グラフにおいて破線で示されたデータは、式(8)のアルゴリズムを使ったときの結果に相当する。
図9Aおよび図9Bを参照すると、図6Aおよび図6Bと比べ、カット数の平均値およびカット数の最大値のいずれも最適解に近づいていることがわかる。ただし、図9Aおよび図9Bの結果でも依然として計算値と最適解の間に差がある。この誤差は、第3ベクトルの要素の値ziを定義する基本式において、変数としてスピンsiではなく、連続変数xiが使われていることに起因していると考えられる。特に、高次の項が増えるほど、ziにおける変数xどうしの積演算は、誤差を大きくする原因となりうる。例えば、1より大きい変数を複数回乗算すると、値が1より著しく大きくなってしまう。
(第2アルゴリズム)
誤差を軽減するために、式(9)の第1アルゴリズムをさらに改良した。具体的には、下記の式(11)のように、z
iにおいて、連続変数x
iに代わって連続変数x
iを符号関数で変換した値sgn(x
i)を代入した。連続変数x
iを符号関数で変換した値sgn(x
i)は、スピンs
iに相当する。
なお、QUBO問題の場合、式(11)のz
i(t+Δt)は、下記の式(12)により表される。
式(11)では、z
iの1階のテンソルを含む項の係数αを定数(例えば、α=1)にしてもよい。式(11)の第2アルゴリズムは、式(8)、(9)とは異なり、ハミルトン方程式を解くものではない。式(11)は、外場によって制御された力学系であると見なすことが可能である。式(11)の第2アルゴリズムでは、高次の目的関数を有するHOMOを扱った場合、z
iのどのスピンどうしの積も-1または1のいずれかの値をとるため、積演算による誤差の発生を防ぐことができる。
上述の式(11)の第2アルゴリズムのように、計算サーバ3が計算するデータは、さらに、si(i=1、2、・・・、N)を要素とする第4ベクトル(s1,s2,・・・,sN)を含んでいてもよい。第1ベクトルのそれぞれの要素を符号関数で変換することにより、第4ベクトルを得ることができる。すなわち、作用演算部51は、イジングモデルのエネルギー式を、すべての項に含まれる変数について偏微分した形式の基本式を使って、第3ベクトルのそれぞれの要素の値を更新するように構成されていてもよい。ここで、基本式の変数として、第1ベクトルの要素または、第1ベクトルの要素を符号関数で変換した第4ベクトルの要素を使うことができる。
図10Aおよび図10Bは、式(11)の第2アルゴリズムを使ってG-setのG22を1000回解いたときの結果を示している。使用されるアルゴリズムの違いを除けば、図10Aおよび図10Bの計算条件(例えば、時間ステップΔt、各係数、xiを-1~1の間に保つために定義されるw)は、図9Aおよび図9Bと同様であるものとする。図10Aは、カット数の平均値を示している。一方、図10Bは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Cmaxの定義は、図6Aおよび図6B、図9Aおよび図9Bと同様である。また、図10Aおよび図10Bの両グラフにおいて実線で示されたデータは、式(11)の第2アルゴリズムの適用時の結果に相当する。一方、図10Aおよび図10Bの両グラフにおいて破線で示されたデータは、式(8)のアルゴリズムを使ったときの結果に相当する。
図10Aおよび図10Bを参照すると、カット数の平均値およびカット数の最大値のいずれにおいても、図9Aおよび図9Bと比べて最適解により近い値が得られることがわかる。図10Bを参照すると、式(11)の第2アルゴリズムを使うことにより、カット数の最大値13359が得られることがわかる。
(第3アルゴリズム)
式(9)の第1アルゴリズムをさらに下記の式(13)のように変形してもよい。
なお、QUBO問題の場合、式(13)のz
i(t+Δt)は、下記の式(14)により表される。
式(13)の第3アルゴリズムは、外力に相当する項fiの計算方法が上述の各例とは異なっている。(13)4番目の式を使って計算された値ziを符号関数で変換し、1で規格化している。すなわち、作用演算部51は、第1ベクトルの要素を変数として計算した基本式の値(zi)を第1関数で変換した値に基づいて第3ベクトルのそれぞれの要素を更新するように構成されていてもよい。第1関数として、例えば、符号関数を使うことができる。ただし、後述するように、その他の関数を第1関数として使ってもよい。
また、式(13)では、係数cの代わりに、関数g(t)が使われている。一般に、第3ベクトルの要素の値z
iの計算結果への寄与度は、問題によって異なる。しかし、式(13)では、第3ベクトルの要素の値z
iが1で規格化されているため、問題ごとに係数cの値を決定する必要がなくなる。関数g(t)として、例えば、下記の式(15)を使うことができる。
式(15)の関数は、更新回数に応じて、単調増加してから単調減少する。ただし、上述の式(15)は一例にすぎず、g(t)として、p(t)をパラメータとする、これとは異なる関数を使ってもよい。すなわち、作用演算部51は、第1係数pをパラメータとする第2関数を乗じることによって第3ベクトルのそれぞれの要素を更新するように構成されていてもよい。
図11Aおよび図11Bは、式(13)の第3アルゴリズムを使ってG-setのG22を1000回解いたときの結果を示している。使用されるアルゴリズムの違いを除けば、図11Aおよび図11Bの計算条件(例えば、時間ステップΔt、使われる係数、xiを-1~1の間に保つために定義されるw)は、図9Aおよび図9Bと同様であるものとする。図11Aは、カット数の平均値を示している。一方、図11Bは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Cmaxの定義は、図6Aおよび図6B、図9Aおよび図9Bと同様である。また、図11Aおよび図11Bの両グラフにおいて実線で示されたデータは、式(13)の第3アルゴリズムの適用時の結果に相当する。一方、図11Aおよび図11Bの両グラフにおいて破線で示されたデータは、式(8)のアルゴリズムを使ったときの結果に相当する。
図11Aおよび図11Bを参照すると、カット数の平均値およびカット数の最大値のいずれにおいても、式(8)のアルゴリズムと比べて最適解により近い値が得られることがわかる。図11Bを参照すると、式(13)の第3アルゴリズムを使うことにより、最大値13359が得られていることがわかる。
(変形例)
式(9)、式(11)および式(13)のアルゴリズムでは、基本式(ziの式)の1階のテンソルを含む項の係数αを定数係数(例えば、α=1)として計算を行ってもよい。また、式(9)、式(11)および式(13)のアルゴリズムでは、基本式(ziの式)の1階のテンソルを含む項の係数αとして更新回数に応じて単調減少または単調増加する係数を使ってもよい。この場合、基本式の1階のテンソルを含む項は、更新回数に応じて単調減少または単調増加する。
上述の式(9)の第1アルゴリズムおよび式(11)の第2アルゴリズムには、係数cが含まれている。係数cをJ
(2)行列の最大固有値の逆数に近い値にしたい場合、J
(2)行列の最大固有値を計算するか、J
(2)行列の最大固有値の見積もりを行う必要がある。最大固有値の計算は、必要な計算量が大きい。一方、最大固有値の見積もりは、値の正確性が保証されていない。そこで、係数cの代わりに上述の式(15)のような、更新回数に応じて値が変動する関数を使うことができる。また、係数cの代わりに、下記の式(16)のように、第1のベクトル(x
1、x
2、・・・、x
N)および第3ベクトル(z
1、z、・・・、z
N)に基づいて計算される、近似値c1を使ってもよい。
式(16)を参照すると、分母と分子がいずれもベクトルのノルムとなっている。式(16)のように、ベクトルのノルムとして、ベクトルの各要素の2乗和の平方根である、L2ノルムを使うことができる。ただし、ベクトルの要素の絶対値の和である、L1ノルムなど、その他の定義によるノルムを使ってもよい。
すなわち、更新部50は、第1ベクトルのノルムを第3ベクトルのノルムで除算した第3係数c1を計算し、第2ベクトルの要素に、第1係数p(t+Δt)と更新された第1ベクトルの対応する要素との積を重み付けした値ならびに、第3ベクトルの対応する要素を第3係数c1で重み付けした値を加算することによって第2ベクトルの要素を更新するように構成されていてもよい。
さらに、係数cの代わりに、下記の式(17)のような内積によって定義される近似値c´1を使ってもよい。
すなわち、更新部50は、第1ベクトルどうしの内積を、第1ベクトルと第3ベクトルの内積の絶対値で除算した第3係数c´1を計算し、第2ベクトルの要素に、第1係数p(t+Δt)と更新された第1ベクトルの対応する要素との積を重み付けした値ならびに、第3ベクトルの対応する要素を第3係数c´1で重み付けした値を加算することによって第2ベクトルの要素を更新するように構成されていてもよい。
近似値c1、c´1は、各計算タイミングにおける第1ベクトル(x1,x2,・・・,xN)および第3ベクトル(z1,z2,・・・,zN)の値に基づいて計算されるため、係数cのような定数ではなく、動的に制御される係数となる。なお、第1ベクトル(x1,x2,・・・,xN)および第3ベクトル(z1,z2,・・・,zN)については、変数の更新処理で計算されたものを利用することができるため、近似値c1、c´1を計算したとしても、計算量が大幅に増えることはない。局所磁場のないイジング問題において、(x1,x2,・・・,xN)がJ(2)の最大固有値に対応する固有ベクトルであるとき、近似値c1、c´1はJ(2)の最大固有値の逆数に等しくなる。また、(x1,x2,・・・,xN)が固有ベクトルからずれていると、近似値c1、c´1はJ(2)の最大固有値の逆数より大きな値となるため、解の収束が早まる。
図12Aおよび図12Bは、式(9)の第1アルゴリズムにおいて、係数cの代わりに近似値c1を使い、G-setのG22を1000回解いたときの結果を示している。図13Aおよび図13Bは、式(11)の第2アルゴリズムにおいて、係数cの代わりに近似値c1を使い、G-setのG22を1000回解いたときの結果を示している。
図12Aおよび図13Aは、カット数の平均値を示している。一方、図12Bおよび図13Bは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Cmaxの定義は、上述の各グラフと同様である。図12A、図12B、図13Aおよび図13Bにおいて破線で示されたデータは、式(8)のアルゴリズムを使ったときの結果を示している。
図12A、図12B、図13Aおよび図13Bを参照すると、カット数の平均値およびカット数の最大値のいずれにおいても、式(8)のアルゴリズムと比べて最適解により近い値が得られることがわかる。特に、式(11)の第2アルゴリズムにおいて、カット数の最大値13359が得られていることがわかる。
なお、式(9)の第1アルゴリズムおよび式(11)の第2アルゴリズムでは、近似値c1およびc´1の代わりに下記の式(18)で定義される近似値c2またはc´2を使ってもよい。
すなわち、更新部50は、第1ベクトルのそれぞれの要素を符号関数で変換した第4ベクトルのノルムを、第3ベクトルのノルムで除算した第3係数c2を計算し、第2ベクトルの要素に、第1係数p(t+Δt)と更新された第1ベクトルの対応する要素との積を重み付けした値ならびに、第3ベクトルの対応する要素を第3係数c2で重み付けした値を加算することによって第2ベクトルの要素を更新するように構成されていてもよい。
また、更新部50は、第1ベクトルのそれぞれの要素を符号関数で変換した第4ベクトルどうしの内積を、第4ベクトルと第3ベクトルの内積の絶対値で除算した第3係数c´2を計算し、第2ベクトルの要素に、第1係数p(t+Δt)と更新された第1ベクトルの対応する要素との積を重み付けした値ならびに、第3ベクトルの対応する要素を第3係数c´2で重み付けした値を加算することによって第2ベクトルの要素を更新するように構成されていてもよい。
式(18)の第3ベクトル(z1,z2,・・・,zN)については、アルゴリズムで計算されたものを使うことができるため、近似値c2、c´2を求めたとしても、計算量が大幅に増えることはない。
アルゴリズムの実行中における各ベクトルの値を使って、近似値c1、c´1、c2、c´2を計算すると、各計算タイミングによって値が激しく変動することがある。近似値c1、c´1、c2、c´2の変動を抑制するため、近似値c1、c´1、c2、c´2の代わりに近似値c1、c´1、c2、c´2を所定の規則に基づいて変換した値を使ってもよい。例えば、所定の規則として下記の式(19)を使うことができる。
ここで、γに1より小さい値を設定することができる。式(19)のc(t+Δt)には、例えば、上述の式(16)~(18)によって計算された近似値が代入される。c(t+Δt)を、各計算タイミングにおいて振動成分を含む信号をサンプリングした値であるとみなすと、d(t+Δt)は、c(t+Δt)が一定帯域のローパスフィルタを通過した後の値に相当するといえる。
すなわち、更新部50は、第3係数(近似値c1、c´1、c2、c´2のいずれか)がローパスフィルタを通過した後の値である第4係数を計算し、第3係数に代わり、第4係数を用いて第2ベクトルの要素を更新するように構成されていてもよい。
上述では、シミュレーテッド分岐アルゴリズムを用いて、イジングモデルの解を求める例について説明した。ただし、シミュレーテッド分岐アルゴリズムによる求解が可能な組合せ最適化問題は、イジング問題に限られない。シミュレーテッド分岐アルゴリズムを用いることによって、一般的な2値変数の組合せ最適化問題を解くことが可能である。例えば、上述の各アルゴリズムは、目的関数の変数が、a(第1値)と、aより大きいb(第2値)のいずれかをとる2値変数である、組合せ最適化問題に適用することが可能である。また、一定の更新回数の後に目的関数の解を求める場合、符号関数の代わりに、値域がaまたはbの2値である関数f(xi)を使ってもよい。この関数f(xi)がとる値は、変数xiの値をしきい値v(a<v<b)と比較した結果に基づいて決まる。例えば、xi<vであるならば、f(xi)=aとなる。また、v<xiであるならば、f(xi)=bとなる。例えば、xi=vである場合、f(xi)=aまたは、f(xi)=bとなる。ここで、しきい値vの値として、例えば、(a+b)/2を使うことができる。上述の関数f(xi)は、第1ベクトルの要素を第4ベクトルの要素に変換する第1関数として使われてもよい。
また、例えば、上述の式(9)の第1アルゴリズム、式(11)の第2アルゴリズムおよび式(13)の第3アルゴリズムを使った場合、更新によって変数xiがaより小さくなったとき、変数xiの値を{v-w-}に変更する。w-は、0以上(v-a)以下の実数である。また、更新によって変数xiがbより大きくなったとき、変数xiの値を{v+w+}に変更する。w+は、0以上(b-v)以下の実数である。
例えば、w-は、0以上(v-a)以下の予め定められた値であってもよい。w+は、0以上(b-v)以下の予め定められた値であってもよい。
また、w-は、0以上(v-a)以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値であってもよい。また、w+は、0以上(b-v)以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値であってもよい。
また、w-およびw+は、複数の第1変数(xi)のしきい値(v)からの偏差の大きさに対する平均を表す指標値(xave)であってもよい。例えば、指標値(xave)は、複数の要素の第1変数(xi)のしきい値(v)からの偏差の二乗平均平方根または絶対値平均である。
また、w-は、指標値(xave)以上、(v-a)以下の乱数により定まる値であってもよい。また、w+は、指標値(xave)以上、(b-v)以下の乱数により定まる値であってもよい。なお、指標値(xave)が(v-a)を超える場合には、w-は、(v-a)とする。また、指標値が(b-v)を超える場合には、w+は、(b-v)とする。
また、w-およびw+は、増加係数であってもよい。増加係数は、更新処理の初期時刻において0であり、初期時刻から終了時刻まで時間経過に従って増加する。また、w-は、増加係数以上、(v-a)以下の乱数により定まる値であってもよい。
また、w+は、増加係数以上、(b-v)以下の乱数により定まる値であってもよい。なお、増加係数が(v-a)を超える場合には、w-は、(v-a)とする。また、増加係数が(b-v)を超える場合には、w+は、(b-v)とする。
以上から、目的関数の変数が第1値(a)と第2値(b)のいずれかをとる離散変数である組合せ最適化問題に対して、上述の式(9)の第1アルゴリズム、式(11)の第2アルゴリズムおよび式(13)の第3アルゴリズムを適用した場合、第1変数xiおよび第2変数yiの更新処理において、更新部50は、次のような処理を行う。
すなわち、更新部50は、第1変数(xi)が第1値(a)より小さい場合、第1変数(xi)を、第1値(a)以上、しきい値(v)以下の値に変更する。さらに、更新部50は、第1変数(xi)が第2値(b)より大きい場合、第1変数(xi)を、予め定められたしきい値(v)以上、第2値(b)以下の値に変更する。
より具体的には、例えば、更新部50は、第1変数(xi)が第1値(a)より小さい場合、第1変数(xi)を、第1値(a)以上しきい値(v)以下の予め定められた値、または、第1値(a)以上しきい値(v)以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値に変更してもよい。さらに、更新部50は、第1変数(xi)が第2値(b)より大きい場合、第1変数(xi)を、しきい値(v)以上第2値(b)以下の予め定められた値、または、しきい値(v)以上第2値(b)以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値に変更してもよい。
また、例えば、更新部50は、第1変数(xi)が第1値(a)より小さい場合、第1変数(xi)を、しきい値(v)から指標値を減算した値に変更してもよい。さらに、更新部50は、第1変数(xi)が第2値(b)より大きい場合、第1変数(xi)を、しきい値(v)に指標値を加算した値に変更してもよい。この場合、指標値は、複数の要素の第1変数(xi)のしきい値(v)からの偏差の大きさに対する平均を表す。例えば、指標値は、複数の要素の第1変数(xi)のしきい値(v)からの偏差の二乗平均平方根または絶対値平均である。
また、例えば、更新部50は、第1変数(xi)が第1値(a)より小さい場合、第1変数(xi)を、第1値(a)以上、しきい値(v)から指標値を減じた値以上の、乱数により定まる値に変更してもよい。さらに、更新部50は、第1変数(xi)が第2値(b)より大きい場合、第1変数(xi)を、しきい値(v)に指標値を加算した値以上、第2値(b)以下の、乱数により定まる値に変更してもよい。
また、更新部50は、第1変数(xi)が第1値(a)より小さい場合、第1変数(xi)を、しきい値(v)から増加係数を減算した値に変更してもよい。さらに、更新部50は、第1変数(xi)が第2値(b)より大きい場合、第1変数(xi)を、しきい値(v)に増加係数を加算した値に変更してもよい。この場合、増加係数は、初期時刻において0であり、初期時刻から終了時刻まで時間経過に従って増加する。
また、更新部50は、第1変数(xi)が第1値(a)より小さい場合、第1変数(xi)を、第1値(a)以上、しきい値(v)に増加係数を加算した値以下における、乱数により定まる値に変更してもよい。さらに、更新部50は、第1変数(xi)が第2値(b)より大きい場合、第1変数(xi)を、しきい値(v)に増加係数を加算した値以上、第2値(b)以下における、乱数により定まる値に変更してもよい。
ここまでは、シンプレクティック・オイラー法によって実装されたシミュレーテッド分岐アルゴリズムの例と、それぞれのアルゴリズムを使って組合せ最適化問題を計算した結果について説明した。以降では、上述のアルゴリズムの実装例について述べる。
(PCクラスタへの実装例)
はじめに、PCクラスタへ上述のアルゴリズムを実装した例について説明する。PCクラスタとは、複数台のコンピュータを接続し、1台のコンピュータでは得られない計算性能を実現するシステムである。例えば、図1に示した情報処理システム100は、複数台の計算サーバ3およびプロセッサを含んでおり、PCクラスタとして利用することが可能である。例えば、PCクラスタにおいては、MPI(Message Passing Interface)を使うことにより、情報処理システム100のような複数の計算サーバ3にメモリが分散して配置されている構成でも並列的な計算を実行することが可能である。例えば、MPIを使って管理サーバ1の制御プログラム14E、各計算サーバ3の計算プログラム34Bおよび制御プログラム34Cを実装することができる。
PCクラスタで利用するプロセッサ数がQである場合、それぞれのプロセッサに、第1ベクトル(x
1,x
2,・・・,x
N)に含まれる変数x
iのうち、L個の変数の計算を行わせることができる。同様に、それぞれのプロセッサに、第2ベクトル(y
1,y
2,・・・,y
N)に含まれる変数y
iのうち、L個の変数の計算を行わせることができる。すなわち、プロセッサ#j(j=1,2,・・・,Q)は、変数{x
m|m=(j-1)L+1,(j-1)L+2,・・・,jL}および{y
m|m=(j-1)L+1,(j-1)L+2,・・・,jL}の計算を行う。また、プロセッサ#jによる{y
m|m=(j-1)L+1,(j-1)L+2,・・・,jL}の計算に必要な下記の式(20)に示されたテンソルJ
(n)は、プロセッサ#jがアクセス可能な記憶領域(例えば、レジスタ、キャッシュ、メモリなど)に保存されるものとする。
ここでは、それぞれのプロセッサが第1ベクトルおよび第2ベクトルの一定数の変数を計算する場合を説明した。ただし、プロセッサによって、計算する第1ベクトルおよび第2ベクトルの変数の数が異なっていてもよい。例えば、計算サーバ3に実装されるプロセッサによって性能差がある場合、プロセッサの性能に応じて計算対象とする変数の数を決めることができる。
すなわち、情報処理装置(例えば、計算サーバ3)は、複数のプロセッサを備えていてもよい。更新部50は、複数のプロセッサを含んでおり、更新部50の複数のプロセッサのそれぞれは、第1ベクトルの一部の要素の値および第2ベクトルの一部の要素の値を更新するように構成されていてもよい。
変数yiの値を更新するためには、第1ベクトル(x1,x2,・・・,xN)または第1ベクトルの各要素を2値変数に変換した第4ベクトル(s1,s2,・・・,sN)のすべての成分の値が必要となる。2値変数への変換は、例えば、符号関数sgn()を使うことによって行うことができる。そこで、Allgather関数を使い、第1ベクトル(x1,x2,・・・,xN)または第4ベクトル(s1,s2,・・・,sN)のすべての成分の値をQ個のプロセッサに共有させることができる。第1ベクトル(x1,x2,・・・,xN)または第4ベクトル(s1,s2,・・・,sN)については、プロセッサ間での値の共有が必要であるものの、第2ベクトル(y1,y2,・・・,yN)およびテンソルJ(n)については、プロセッサ間での値の共有を行うことは必須ではない。プロセッサ間でのデータの共有は、例えば、プロセッサ間通信を使ったり、共有メモリにデータを保存したりすることによって実現することができる。
プロセッサ#jは、変数{zm|m=(j-1)L+1,(j-1)L+2,・・・,jL}の値を計算する。そして、プロセッサ#jは、計算した{zm|m=(j-1)L+1,(j-1)L+2,・・・,jL}の値に基づき、変数{ym|m=(j-1)L+1,(j-1)L+2,・・・,jL}を更新する。
上述の各式に示したように、ベクトル(z1,z2,・・・,zN)の計算では、テンソルJ(n)と、ベクトル(x1,x2,・・・,xN)または(s1,s2,・・・,sN)との積の計算を含む、積和演算が必要である。積和演算は、上述のアルゴリズムにおいて最も計算量の大きい処理であり、計算速度の向上においてボトルネックとなりうる。そこで、PCクラスタの実装では、積和演算を、Q=N/L個のプロセッサに分散して並列的に実行し、計算時間の短縮をはかることができる。
すなわち、情報処理装置(例えば、計算サーバ3)は、複数のプロセッサを備えていてもよい。作用演算部51は、複数のプロセッサを含み、作用演算部51の複数のプロセッサのそれぞれは、第3ベクトルの一部の要素を更新するように構成されていてもよい。更新部50は、複数のプロセッサを含み、更新部50の複数のプロセッサのそれぞれは、第1ベクトルの一部の要素および第2ベクトルの一部の要素を更新するように構成されていてもよい。
図14は、マルチプロセッサ構成の例を概略的に示している。図14の複数の計算ノードは、例えば、情報処理システム100の複数の計算サーバ3に相当する。また、図14の高速リンクは、例えば、情報処理システム100のケーブル4a~4cおよびスイッチ5によって形成された計算サーバ3間のインターコネクトに相当する。図14の共有メモリは、例えば、共有メモリ32に相当する。図14のプロセッサは、例えば、各計算サーバ3のプロセッサ33A~33Dに相当している。なお、図14には複数の計算ノードが示されているが、単一計算ノードの構成を用いることを妨げるものではない。
図14には、各構成要素に配置されるデータおよび構成要素間で転送されるデータが示されている。各プロセッサでは、変数xi、(si)、yi、ziの値が計算される。また、プロセッサと共有メモリ間では、変数xiまたはsiが転送される。各計算ノードの共有メモリには、例えば、第1ベクトル(x1,x2,・・・,xN)、第2ベクトル(y1,y2,・・・,yN)のL個の変数、およびテンソルJ(n)の一部が保存される。なお、式(11)のアルゴリズムを実行する場合、各計算ノードの共有メモリには、第1ベクトル(x1,x2,・・・,xN)に代わって第4ベクトル(s1,s2,・・・,sN)が保存されていてもよい。そして、計算ノード間を接続する高速リンクでは、例えば、第1ベクトル(x1,x2,・・・,xN)が転送される。Allgather関数を使った場合、各プロセッサで変数yiおよびziを更新するために、第1ベクトル(x1,x2,・・・,xN)の全要素が必要だからである。なお、式(11)のアルゴリズムにしたがって変数ziを更新する場合には、各プロセッサは第4ベクトル(s1,s2,・・・,sN)の全要素にアクセスする必要がある。このため、高速リンクでは、第4ベクトル(s1,s2,・・・,sN)が転送されてもよい。
ただし、図14に示したデータの配置および転送は一例にしかすぎない。例えば、各プロセッサが積和演算を含む{zm|m=(j-1)L+1,(j-1)L+2,・・・,jL}の計算を並列的に実行しているのであれば、それぞれのプロセッサと共有メモリ間および計算ノード間で変数ziの値を転送し、共有されたベクトル(z1,z2,・・・,zN)を参照し、変数yiの値を計算してもよい。このように、PCクラスタにおけるデータの配置方法、転送方法および並列化の実現方法については、特に問わない。
すなわち、情報処理装置(例えば、計算サーバ3)は、複数のプロセッサからアクセス可能に構成されている共有メモリを備えていてもよい。この場合、更新部50は、更新された後の第1ベクトルの要素または、更新された後の第1ベクトルのそれぞれの要素を2値変数に変換した第4ベクトルを共有メモリに保存することができる。
次に、PCクラスタに上述の各アルゴリズムを実行させたときにおける結果について説明する。図15Aおよび図15Bは、PCクラスタを使ってN=3600の(局所磁場のない)全結合イジング問題を解いたときの結果を示している。図15Aおよび図15Bの全結合イジング問題では、結合係数の行列Jの各成分の値は、[-1,1]の範囲の一様乱数に設定した。また、合計時間ステップ数は、10000とした。図15Aは、各アルゴリズムを使って10回全結合イジング問題を解いたときにおけるエネルギーEIsingの平均値を示している。図15Bは、各アルゴリズムを使って10回全結合イジング問題を解いたときにおける計算時間の平均値を秒単位で示している。
図15Aおよび図15Bの棒グラフには、左側から右側に向かって、
(i) 式(8)のアルゴリズムが使われた場合、
(ii) 式(9)の第1アルゴリズムが使われた場合、
(iii) 式(11)の第2アルゴリズムが使われた場合、
(iv) 式(13)の第3アルゴリズムが使われた場合、
(v) 式(9)の第1アルゴリズムにおいて、係数cの代わりに近似値c1が使われた場合、
(vi) 式(11)の第2アルゴリズムにおいて、係数cの代わりに近似値c1が使われた場合、
の6のケースの結果が示されている。
また、図15Aおよび図15Bの棒グラフには、各ケースのそれぞれについて、左側から右側に向かって、プロセッサ数Qが1および36の場合における結果が示されている。式(8)のアルゴリズムが使われた場合、式(9)の第1アルゴリズムが使われた場合、式(11)の第2アルゴリズムが使われた場合において、係数cとして定数0.5D√(3/N)を用いた。
図15Aの棒グラフを参照すると、(ii)~(vi)のケースでは、(i)のケースと比べ、エネルギーEIsingの平均値が低くなっており、最適解により近い解が得られやすいことがわかる。また、図15Bの棒グラフを参照すると、マルチプロセッサ構成による並列計算を行うことにより、計算時間が著しく短縮されることがわかる。
(GPUへの実装例)
GPU(Graphics Processing Unit)を使って上述の各アルゴリズムの計算を行ってもよい。図16は、GPUを使った構成の例を概略的に示している。図16には、互いに高速リンクで接続された複数のGPUが示されている。それぞれのGPUには、共有メモリにアクセス可能な複数のコアが搭載されている。また、図16の構成例では、複数のGPUが高速リンクを介して接続されており、GPUクラスタを形成している。例えば、GPUが図1のそれぞれの計算サーバ3に搭載されている場合、高速リンクは、ケーブル4(4a~4c)およびスイッチ5によって形成された計算サーバ3間のインターコネクトに相当する。なお、図16の構成例では、複数のGPUが使われているが、ひとつのGPUを使った場合にも、並列的な計算を実行することが可能である。すなわち、図16のそれぞれのGPUは、図14のそれぞれの計算ノードに相当する計算を実行できる。すなわち、情報処理装置(計算サーバ3)のプロセッサは、Graphics Processing Unit(GPU)のコアであってもよい。
GPUにおいて、変数xiおよびyi、ならびにテンソルJ(n)はデバイス変数として定義される。GPUは、変数yiの更新に必要なテンソルJ(n)と第1ベクトル(x1,x2,・・・,xN)または第4ベクトル(s1,s2,・・・,sN)の積を、行列ベクトル積関数によって並列的に計算することができる。なお、行列とベクトルの積演算を繰り返し実行することにより、テンソルとベクトルの積を求めることができる。また、第1ベクトル(x1,x2,・・・,xN)の計算と、第2ベクトル(y1,y2,・・・,yN)のうち、積和演算以外の部分については、それぞれのスレッドにi番目の要素(xi,yi)の更新処理を実行させ、処理の並列化を実現することができる。
図17Aおよび図17Bは、GPUを使ってN=3600の全結合イジング問題を解いたときの結果を示している。図17Aの棒グラフは、各アルゴリズムを使って10回全結合イジング問題を解いたときにおけるエネルギーEIsingの平均値を示している。図17Bの棒グラフは、各アルゴリズムを使って10回全結合イジング問題を解いたときにおける計算時間の平均値を秒単位で示している。
図17Aおよび図17Bの棒グラフには、左側から右側に向かって、式(8)のアルゴリズムが使われた場合、式(9)の第1アルゴリズムが使われた場合、式(11)の第2アルゴリズムが使われた場合、式(13)の第3アルゴリズムが使われた場合の4つのケースにおける結果が示されている。式(8)のアルゴリズムが使われた場合、式(9)の第1アルゴリズムが使われた場合、式(11)の第2アルゴリズムが使われた場合において、係数cとして定数0.5D√(3/N)を用いた。いずれの結果においても合計時間ステップ数は、10000である。また、それぞれのアルゴリズムについて、左側が1計算ノードのPCクラスタ、右側が1GPUを用いたときの結果が示されている。
図17Aの棒グラフを参照すると、式(9)、式(11)および式(13)のアルゴリズムが使われた場合には、式(8)のアルゴリズムが使われた場合と比べ、エネルギーEIsingの平均値が低くなっており、最適解により近い解が得られやすいことがわかる。また、図17Bの棒グラフを参照すると、GPUによる並列計算を行うことにより、1計算ノードのPCクラスタと比べて計算時間を著しく短縮できることがわかる。これは、GPUの計算の並列度が、一般的なCPUと比べて高いことに起因する。
(第1アルゴリズムの処理フローの第1例)
図18は、第1アルゴリズムを実行する場合における情報処理システム100の処理の流れの第1例を示す。情報処理システム100は、式(9)および式(10)に示す第1アルゴリズムを用いて最適化問題を解く場合、例えば、図18に示す流れで処理を実行する。
まず、S101において、更新部50は、パラメータを設定する。具体的には、更新部50は、N×N個の結合係数を含む行列であるJ、および、N個の局所磁場を表す局所磁場係数を含む配列であるhを設定する。なお、更新部50は、HOBO問題を解く場合には、Jおよびhに代えて、Nn個の作用係数を含むn次テンソルであるJ(n)を設定する。この場合、nは、HOBO問題の目的関数の変数の次数を表す。さらに、更新部50は、係数であるD、係数であるc、単位時間を表すΔt、終了時刻を表すT、関数であるp(t)、および、関数であるα(t)を設定する。p(t)およびα(t)は、t=初期時刻(例えば0)で0、t=終了時刻(T)で1となる増加関数である。更新部50は、J、hをユーザからの受け取った情報に応じて設定する。更新部50は、D、c、Δt、T、p(t)およびα(t)を、ユーザから受け取ったパラメータに応じて設定してもよいし、予め決定されており変更できないパラメータを設定してもよい。
続いて、S102において、更新部50は、変数を初期化する。具体的には、更新部50は、時刻を表す変数であるtを初期時刻(例えば、0)に初期化する。さらに、更新部50は、N個の第1変数(x1(t)~xN(t))のそれぞれおよびN個の第2変数(y1~yN)のそれぞれに、ユーザから受け取った初期値、予め定められた固定値、または、乱数を代入する。
続いて、更新部50は、S103とS118との間のループ処理を、tがTより大きくなるまで繰り返す。1回のループ処理において、更新部50は、対象時刻(t+Δt)におけるN個の第1変数(x1(t+Δt)~xN(t+Δt))を、直前時刻(t)におけるN個の第2変数(y1(t)~yN(t))に基づき算出する。また、1回のループ処理において、更新部50は、対象時刻(t+Δt)におけるN個の第2変数(y1(t+Δt)~yN(t+Δt))を、直前時刻(t)におけるN個の第1変数(x1(t)~xN(t))に基づき算出する。
なお、直前時刻(t)は、対象時刻(t+Δt)より単位時間(Δt)前の時刻である。すなわち、更新部50は、S103とS118との間のループ処理を繰り返すことにより、N個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を、初期時刻(t=0)から終了時刻(t=T)まで単位時間(Δt)毎に順次に更新する。
S104において、更新部50は、直前時刻(t)におけるN個の第1変数(x
1(t)~x
N(t))の大きさの平均を表す指標値(x
ave)を算出する。例えば、指標値(x
ave)は、直前時刻におけるN個の第1変数(x
1(t)~x
N(t))の二乗平均平方根または絶対値平均である。例えば、更新部50は、二乗平均平方根を算出する場合、式(21-1)に示す演算を実行する。また、例えば、更新部50は、絶対値平均を算出する場合、式(21-2)に示す演算を実行する。なお、更新部50は、後述するS109において、指標値(x
ave)を用いない場合、S104の処理を実行しない。
続いて、更新部50は、S105とS111との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。iは、1からNまでの整数であり、N個の要素のうちの処理対象を表すインデックスである。N個の要素のそれぞれは、第1変数(xi(t))および第2変数(yi(t))が対応付けられる。S105とS111との間のループ処理において、更新部50は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S106において、更新部50は、対象要素の対象時刻(t+Δt)における第1変数(x
i(t+Δt))を、対象要素の直前時刻(t)における第1変数(x
i(t))に、対象要素の直前時刻(t)における第2変数(y
i(t))と予め定められた定数(D)と単位時間(Δt)とを乗算した値を加算することにより算出する。具体的には、更新部50は、式(22)を算出する。
続いて、S107において、更新部50は、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が、予め定められた第2値(+1)より大きいか否かを判断する。本例において、第2値は、+1である。第2値は、連続量である第1変数(xi(t))の単位量である。更新部50は、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が第2値以下である場合、処理をS111に進める(S107のNo)。更新部50は、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が第2値より大きい場合、処理をS108に進める。
続いて、S108において、更新部50は、0以上、第2値(例えば+1)以下の範囲内の所定区間で一様な確率で発生する乱数(r)を生成する。なお、更新部50は、後述するS109において、乱数(r)を用いない場合、S108の処理を実行しない。
続いて、S109において、更新部50は、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))の制約処理をする。具体的には、更新部50は、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、符号を変更せずに、絶対値を0以上、第2値以下の値に変更する。本例においては、更新部50は、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))が第1値(a)である-1より小さい場合、第1変数(xi(t+Δt))を、-1以上、しきい値(v)である0以下の値に変更する。さらに、更新部50は、第1変数(xi(t+Δt))が第2値(b)である+1より大きい場合、第1変数(xi(t+Δt))を、0以上、+1以下の値に変更する。
例えば、更新部50は、第1変数(xi(t+Δt))が-1より小さい場合、第1変数(xi(t+Δt))を、-1以上0以下の予め定められた値、または、S108で生成した乱数(r)に変更してもよい。さらに、更新部50は、第1変数(xi(t+Δt))が+1より大きい場合、第1変数(xi(t+Δt))を、0以上+1以下の予め定められた値、または、0からS108で生成した乱数(r)を減じた値に変更してもよい。
また、例えば、更新部50は、第1変数(xi(t+Δt))が-1より小さい場合、第1変数(xi(t+Δt))を、0からS104で算出した指標値(xave)を減算した値に変更してもよい。さらに、更新部50は、第1変数(xi(t+Δt))が+1より大きい場合、第1変数(xi(t+Δt))を、0に指標値(xave)を加算した値に変更してもよい。
また、例えば、更新部50は、第1変数(x
i(t+Δt))が-1より小さい場合、第1変数(x
i(t+Δt))を、-1以上、0から指標値(x
ave)を減じた値以上の、S108で算出した乱数(r)により定まる値に変更してもよい。さらに、更新部50は、第1変数(x
i(t+Δt))が+1より大きい場合、第1変数(x
i(t+Δt))を、0に指標値(x
ave)を加算した値以上、+1以下の、S108で算出した乱数(r)により定まる値に変更してもよい。この場合、更新部は、例えば、式(23)を算出する。
また、更新部50は、第1変数(xi(t+Δt))が-1より小さい場合、第1変数(xi(t+Δt))を、0から増加係数を減算した値に変更してもよい。増加関数は、0から+1以下の初期時刻(t=0)から終了時刻(t=T)まで時間経過に従って増加する増加係数に変更してもよい。この場合、更新部50は、初期時刻(t=0)で0、終了時刻(t=T)で第2値(例えば+1)となる、時刻(t)を変数とする一次関数、または、この一次関数の平方根により増加係数を算出する。さらに、更新部50は、第1変数(xi(t+Δt))が+1より大きい場合、第1変数(xi(t+Δt))を、0に増加係数を加算した値に変更してもよい。
また、更新部50は、第1変数(xi(t+Δt))が-1より小さい場合、第1変数(xi(t+Δt))を、-1以上、0に増加係数を加算した値以下における、S108で精鋭した乱数(r)により定まる値に変更してもよい。さらに、更新部50は、第1変数(xi(t+Δt))が+1より大きい場合、第1変数(xi(t+Δt))を、0に増加係数を加算した値以上、+1以下における、S108で生成した乱数(r)により定まる値に変更してもよい。この場合、更新部50は、式(23)に含まれる指標値(xave)を増加係数に代えた式により、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を算出する。
続いて、S128において、更新部50は、対象要素の直前時刻(t)における第2変数(yi(t))の制約処理をする。具体的には、更新部50は、対象要素の直前時刻(t)における第2変数(yi(t))を、0、予め定められた値、または、乱数に応じた値に変更する。更新部50は、乱数に応じた値に変更する場合、例えば、更新部50は、対象要素の直前時刻(t)における第2変数(yi(t))を、予め定められた範囲(例えば、-0.1以上、+0.1以下)内で一様な確率で発生する乱数に変更する。更新部50は、S110を終えると、処理をS111に進める。
更新部50は、以上のようなS105とS111との間のループ処理をN回実行することにより、次の処理を実行する。すなわち、更新部50は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、対象要素の直前時刻(t)における第2変数(yi(t))に基づき更新する。さらに、更新部50は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が、第2値より大きい場合、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、符号を変更せずに、絶対値を0以上、第2値以下に変更する。さらに、更新部50は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が、第2値より大きい場合、対象要素の直前時刻(t)における第2変数(yi(t))を、0、予め定められた値または乱数に応じた値に変更する。
更新部50は、S105とS111との間のループ処理をN回実行した場合、処理をS112に進める。
続いて、更新部50は、S112とS116との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。S112とS116との間のループ処理において、更新部50は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S113において、更新部50は、N個の要素のそれぞれの対象時刻(t+Δt)における第1変数(x1(t+Δt)~xN(t+Δt))と、対象要素とN個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値(zi(t+Δt))を算出する。QUBO問題の場合、作用係数は、Jに含まれる結合係数およびhに含まれる局所磁場係数である。HOBO問題の場合、作用係数は、J(n)に含まれる。
QUBO問題の場合、更新部50は、式(24)を算出する。
また、HOBO問題の場合、更新部50は、式(25)を算出する。
続いて、S114において、更新部50は、対象時刻(t+Δt)における更新値(z
i(t+Δt))に係数(c)を乗算することにより、外力(f
i(t+Δt))を算出する。具体的には、更新部50は、式(26)を算出する。
続いて、S115において、更新部50は、対象要素の対象時刻(t+Δt)における第2変数(y
i(t+Δt))を、対象要素の直前時刻(t)における第2変数(y
i(t))に、外力(f
i(t+Δt))と対象要素の対象時刻(t+Δt)における第1変数(x
i(t+Δt))とに基づく値に単位時間(Δt)を乗算した値を、加算することにより、算出する。具体的には、更新部50は、式(27)を算出する。
更新部50は、以上のようなS112とS116との間のループ処理をN回実行することにより、次の処理を実行する。すなわち、更新部50は、N個の要素のそれぞれについて、対象時刻(t+Δt)における第2変数(yi(t+Δt))を、N個の要素のそれぞれの対象時刻(t+Δt)における第1変数(x1(t+Δt)~xN(t+Δt))に基づき更新する。
更新部50は、S112とS116との間のループ処理をN回実行した場合、処理をS117に進める。S117において、更新部50は、直前時刻(t)に単位時間(Δt)を加算して、対象時刻(t+Δt)を更新する。S118において、更新部50は、S104からS117までの処理を、tが終了時刻(T)を超えるまで繰り返す。そして、更新部50は、tが終了時刻(T)より大きくなった場合、本フローを終了する。
そして、更新部50は、N個の要素のそれぞれについて、終了時刻(t=T)における第1変数(xi(T))の符号に応じて、対応するスピンの値を算出する。例えば、更新部50は、終了時刻(t=T)における第1変数(xi(T))の符号が負である場合、対応するスピンを-1とし、正である場合、対応するスピンを+1とする。そして、更新部50は、算出した複数のスピンの値を組合せ最適化問題の解として出力する。
図18に示すフローチャートに従った処理を実行することにより、更新部50は、N個の要素のそれぞれについて、初期時刻(t=0)から終了時刻(t=T)まで単位時間(Δt)毎の第1変数(xi(t+Δt))および第2変数(yi(t+Δt))を、単位時間毎に順次に、且つ、第1変数(xi(t+Δt))および第2変数(yi(t+Δt))を交互に更新する。また、図18に示すフローチャートに従った処理を実行することにより、更新部50は、単位時間毎に、対象時刻(t+Δt)における第1変数(xi(t+Δt))を算出した後に、対象時刻(t+Δt)における第2変数(yi(t+Δt))を算出する。これにより、更新部50は、シンプレクティック・オイラー法を用いて第1アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
(第1アルゴリズムの処理フローの第2例)
図19は、第1アルゴリズムを実行する場合における情報処理システム100の処理の流れの第2例を示す。なお、図19のフローチャートは、図18で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。
情報処理システム100は、式(9)および式(10)に示す第1アルゴリズムを用いて最適化問題を解く場合、例えば、図19に示す流れで処理を実行してもよい。
まず、S101およびS102において、更新部50は、図18に示す第1例と同様の処理を実行する。続いて、更新部50は、S103とS118との間のループ処理を、図18に示す第1例と同様に実行する。
S104において、更新部50は、図18に示す第1例と同様の処理を実行する。なお、更新部50は、S104の処理をS126の直前で実行してもよい。
続いて、更新部50は、S121とS125との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。S121とS125との間のループ処理において、更新部50は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S122において、更新部50は、N個の要素のそれぞれの直前時刻(t)における第1変数(x1(t)~xN(t))と、対象要素とN個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値(zi(t))を算出する。
QUBO問題の場合、更新部50は、式(28)を算出する。
また、HOBO問題の場合、更新部50は、式(29)を算出する。
続いて、S123において、更新部50は、直前時刻(t)における更新値(z
i(t))に係数(c)を乗算することにより、直前時刻(t)における外力(f
i(t))を算出する。具体的には、更新部50は、式(30)を算出する。
続いて、S124において、更新部50は、対象要素の対象時刻(t+Δt)における第2変数(y
i(t+Δt))を、対象要素の直前時刻(t)における第2変数(y
i(t))に、外力(f
i(t))と対象要素の直前時刻(t)における第1変数(x
i(t))とに基づく値に単位時間(Δt)を乗算した値を、加算することにより、算出する。具体的には、更新部50は、式(31)を算出する。
更新部50は、以上のようなS121とS125との間のループ処理をN回実行することにより、次の処理を実行する。すなわち、更新部50は、N個の要素のそれぞれについて、対象時刻(t+Δt)における第2変数(yi(t+Δt))を、N個の要素のそれぞれの直前時刻(t)における第1変数(x1(t)~xN(t))に基づき更新する。
更新部50は、S121とS125との間のループ処理をN回実行した場合、処理をS126に進める。
続いて、更新部50は、S126とS129との間のループ処理を、i=1からi=Nまでiを1ずつインクリメントしながら繰り返す。S126とS129との間のループ処理において、更新部50は、N個の要素のうちのi番目の要素を、対象要素として処理を実行する。
S127において、更新部50は、対象要素の対象時刻(t+Δt)における第1変数(x
i(t+Δt))を、対象要素の直前時刻(t)における第1変数(x
i(t))に、対象要素の対象時刻(t+Δt)における第2変数(y
i(t+Δt))と予め定められた定数(D)と単位時間(Δt)とを乗算した値を加算することにより算出する。具体的には、更新部50は、式(32)を算出する。
続いて、S107において、更新部50は、図18に示す第1例と同様の処理を実行する。更新部50は、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が第2値以下である場合、処理をS129に進める(S107のNo)。更新部50は、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が第2値より大きい場合、処理をS108に進める。
続いて、S108およびS109において、更新部50は、図18に示す第1例と同様の処理を実行する。
続いて、S110において、更新部50は、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δt))の制約処理をする。具体的には、更新部50は、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δt))を、0、予め定められた値、または、乱数に応じた値に変更する。更新部50は、乱数に応じた値に変更する場合、例えば、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δt))を、予め定められた範囲(例えば、-0.1以上、+0.1以下)内で一様な確率で発生する乱数に変更する。更新部50は、S128を終えると、処理をS129に進める。
更新部50は、以上のようなS126とS129との間のループ処理をN回実行することにより、次の処理を実行する。すなわち、更新部50は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δt))に基づき更新する。さらに、更新部50は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が、第2値より大きい場合、対象要素の対象時刻(t+Δt)における第1変数(xi(t+Δt))を、符号を変更せずに、絶対値を0以上、第2値以下に変更する。さらに、更新部50は、N個の要素のそれぞれについて、対象要素の対象時刻(t+Δt)における第1変数の絶対値(|xi(t+Δt)|)が、第2値より大きい場合、対象要素の対象時刻(t+Δt)における第2変数(yi(t+Δt))を、0、予め定められた値または乱数に応じた値に変更する。
更新部50は、S126とS129との間のループ処理をN回実行した場合、処理をS117に進める。S117において、更新部50は、図18に示す第1例と同様の処理を実行する。S118において、更新部50は、S104からS117までの処理を、tが終了時刻(T)を超えるまで繰り返す。そして、更新部50は、tが終了時刻(T)より大きくなった場合、本フローを終了する。さらに、更新部50は、N個の要素のそれぞれについて、終了時刻(t=T)における第1変数(xi(T))の符号に応じて、対応するスピンの値を算出し、算出した複数のスピンの値を組合せ最適化問題の解として出力する。
図19に示すフローチャートに従った処理を実行することにより、更新部50は、N個の要素のそれぞれについて、初期時刻(t=0)から終了時刻(t=T)まで単位時間(Δt)毎の第1変数(xi(t+Δt))および第2変数(yi(t+Δt))を、単位時間毎に順次に、且つ、第1変数(xi(t+Δt))および第2変数(yi(t+Δt))を交互に更新する。また、図19に示すフローチャートに従った処理を実行することにより、更新部50は、単位時間毎に、対象時刻(t+Δt)における第2変数(yi(t+Δt))を算出した後に、対象時刻(t+Δt)における第1変数(xi(t+Δt))を算出する。これにより、更新部50は、シンプレクティック・オイラー法を用いて第1アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
(第2アルゴリズムの処理フローの第1例)
図20は、第2アルゴリズムを実行する場合における情報処理システム100の処理の流れの第1例を示す。
情報処理システム100は、式(11)および式(12)に示す第2アルゴリズムを用いて最適化問題を解く場合、例えば、図20に示す流れで処理を実行する。なお、図20のフローチャートは、図18で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図18で示した流れで第1アルゴリズムを実行する場合との相違点を説明する。
図20の処理は、S111とS112との間において、S201の処理が追加されている。S201において、更新部50は、N個の要素のそれぞれについて、対象時刻(t+Δt)における第1変数(x
j(t+Δt))の符号(s
j(t+Δt))を算出する。具体的には、j=1~J=Nまでのそれぞれについて、更新部50は、式(33)を算出する。
また、更新部50は、S113に代えて、S202を実行する。S202において、更新部50は、N個の要素のそれぞれの対象時刻(t+Δt)における第1変数の符号(s1(t+Δt)~sN(t+Δt))と、対象要素とN個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値(zi(t+Δt))を算出する。
QUBO問題の場合、更新部50は、式(34)を算出する。
また、HOBO問題の場合、更新部50は、式(35)を算出する。
図20に示すフローチャートに従った処理を実行することにより、更新部50は、シンプレクティック・オイラー法を用いて第2アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
(第2アルゴリズムの処理フローの第2例)
図21は、第2アルゴリズムを実行する場合における情報処理システム100の処理の流れの第2例を示す。
情報処理システム100は、式(11)および式(12)に示す第2アルゴリズムを用いて最適化問題を解く場合、例えば、図21に示す流れで処理を実行することもできる。なお、図21のフローチャートは、図19で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図19で示した流れで第1アルゴリズムを実行する場合との相違点を説明する。
図21の処理は、S121の前において、S211の処理が追加されている。なお、S211は、S104の前に追加されてもよい。
S211において、更新部50は、N個の要素のそれぞれについて、直前時刻(t)における第1変数(x
j(t))の符号(s
j(t))を算出する。更新部50は、式(36)を演算することにより、N個の要素のそれぞれについて、直前時刻(t)における第1変数(x
j(t))の符号(s
j(t))を算出する。具体的には、j=1~J=Nまでのそれぞれについて、更新部50は、式(36)を算出する。
また、更新部50は、S122に代えて、S212を実行する。S212において、更新部50は、N個の要素のそれぞれの直前時刻(t)における第1変数の符号(s1(t)~sN(t))と、対象要素とN個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値(zi(t))を算出する。
QUBO問題の場合、更新部50は、式(37)を算出する。
また、HOBO問題の場合、更新部50は、式(38)を算出する。
図21に示すフローチャートに従った処理を実行することにより、更新部50は、シンプレクティック・オイラー法を用いて第2アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
(第3アルゴリズムの処理フローの第1例)
図22は、第3アルゴリズムを実行する場合における情報処理システム100の処理の流れの第1例を示す。
情報処理システム100は、式(13)および式(14)に示す第3アルゴリズムを用いて最適化問題を解く場合、例えば、図22に示す流れで処理を実行する。なお、図22のフローチャートは、図18で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図18で示した流れで第1アルゴリズムを実行する場合との相違点を説明する。
更新部50は、S114に代えて、S301を実行する。S301において、更新部50は、対象時刻(t+Δt)における更新値(z
i(t+Δt))の符号に、予め定められた関数により定まる係数を乗算することにより、対象時刻(t+Δt)における外力(f
i(t+Δt))を算出する。具体的には、更新部50は、式(39)を算出する。
なお、g(t+Δt)は、式(40)により表される。
図22に示すフローチャートに従った処理を実行することにより、更新部50は、シンプレクティック・オイラー法を用いて第3アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
(第3アルゴリズムの処理フローの第2例)
図23は、第3アルゴリズムを実行する場合における情報処理システム100の処理の流れの第2例を示す。
情報処理システム100は、式(13)および式(14)に示す第3アルゴリズムを用いて最適化問題を解く場合、例えば、図23に示す流れで処理を実行することもできる。なお、図23のフローチャートは、図19で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図19で示した流れで第1アルゴリズムを実行する場合との相違点を説明する。
更新部50は、S123に代えて、S311を実行する。S311において、更新部50は、直前時刻(t)における更新値(z
i(t))の符号に、予め定められた関数により定まる係数を乗算することにより、直前時刻(t)における外力(f
i(t))を算出する。具体的には、更新部50は、式(41)を演算することにより、直前時刻(t)における外力(f
i(t))を算出する。
図23に示すフローチャートに従った処理を実行することにより、更新部50は、シンプレクティック・オイラー法を用いて第3アルゴリズムに従った演算を実行して、終了時刻(t=T)におけるN個の第1変数(x1(t)~xN(t))およびN個の第2変数(y1(t)~yN(t))を算出することができる。
(半導体装置への回路化)
図24は、式(9)の第1アルゴリズム、式(11)の第2アルゴリズムまたは式(13)の第3アルゴリズムを実行することにより、組合せ最適化問題の解を算出する演算装置60の構成を示す図である。演算装置60は、FPGA、ゲートアレイまたは特定用途向け集積回路等であって、半導体装置に回路化して実現される。
演算装置60は、演算回路61と、入力回路62と、出力回路63と、設定回路64とを備える。
演算回路61は、時刻を表すパラメータであるtを、初期時刻(例えば0)から終了時刻まで、単位時間(Δt)ずつ増加させる。演算回路61は、N個の要素(仮想的な粒子)のそれぞれに対応付けられた、第1変数(xi)および第2変数(yi)を算出する。第1変数(xi)は、対応する要素(仮想的な粒子)の位置を表す。第2変数(yi)は、対応する要素(仮想的な粒子)の運動量を表す。
演算回路61は、初期時刻から終了時刻まで単位時間毎に、N個の第1変数(xi)のそれぞれおよびN個の第2変数(yi)のそれぞれを、単位時間毎に順次に且つ第1変数(xi)と第2変数(yi)とを交互に算出する。より具体的には、演算回路61は、初期時刻から終了時刻まで単位時間毎に、式(9)、式(11)または式(13)のアルゴリズムで示された演算を実行する。そして、演算回路61は、終了時刻におけるN個の第1変数(xi)のそれぞれの値(すなわち、N個の仮想的な粒子のそれぞれの位置)を2値化することにより、組合せ最適化問題の解を算出する。
入力回路62は、演算回路61による演算処理に先だって、初期時刻におけるN個の第1変数(xi)およびN個の第2変数(yi)のそれぞれの初期値(すなわち、複数の仮想的な粒子のそれぞれの初期位置および初期運動量)を取得して、演算回路61に与える。出力回路63は、演算回路61による演算処理の終了後に、組合せ最適化問題の解を演算回路61から取得する。そして、出力回路63は、取得した解を出力する。設定回路64は、演算回路61による演算処理に先だって、各パラメータを演算回路61に対して設定する。
(第1アルゴリズムの処理を実行する演算回路61)
図25は、式(9)の第1アルゴリズムを実行する演算回路61のブロック構成を示す図である。演算回路61は、式(9)の第1アルゴリズムを実行する場合、図25に示すように構成される。なお、図25の説明において、直前時刻は、t1と表される。また、対象時刻は、t2と表される。図26および図27でも同様である。
式(9)の第1アルゴリズムを実行する演算回路61は、Xメモリ66と、Yメモリ67と、作用演算回路68と、更新回路69と、制御回路70とを備える。
Xメモリ66は、直前時刻(t1)におけるN個の第1変数(xi(t1))を記憶する。Xメモリ66は、時刻の更新に伴って、記憶している直前時刻(t1)におけるN個の第1変数(xi(t1))が書き換えられる。すなわち、Xメモリ66は、対象時刻(t2)におけるN個の第1変数(xi(t2))が算出された場合、算出された対象時刻(t2)におけるN個の第1変数(xi(t2))が、新たな直前時刻(t1)におけるN個の第1変数(xi(t1))として書き込まれる。設定回路64は、演算に先だって、初期時刻におけるN個の第1変数xiをXメモリ66に書き込む。
Yメモリ67は、直前時刻(t1)におけるN個の第2変数(yi(t1))を記憶する。Yメモリ67は、時刻の更新に伴って、記憶している直前時刻(t1)におけるN個の第2変数(yi(t1))が書き換えられる。すなわち、Yメモリ67は、対象時刻(t2)におけるN個の第2変数(yi(t2))が算出された場合、算出された対象時刻(t2)におけるN個の第2変数(yi(t2))が、新たな直前時刻(t1)におけるN個の第2変数(yi(t1))として書き込まれる。設定回路64は、演算に先だって、初期時刻におけるN個の第2変数yiをYメモリ67に書き込む。
作用演算回路68は、直前時刻(t1)におけるN個の第1変数(xj(t1))をXメモリ66から取得する。そして、作用演算回路68は、N個の要素のそれぞれについて、直前時刻(t1)における更新値(zi(t1))を算出する。
更新回路69は、N個の要素のそれぞれについて、直前時刻(t1)における更新値(zi(t1))を作用演算回路68から取得する。さらに、更新回路69は、N個の要素のそれぞれについて、直前時刻(t1)における第1変数(xi(t1))をXメモリ66から取得し、直前時刻(t1)における第2変数(yi(t1))をYメモリ67から取得する。そして、更新回路69は、N個の要素のそれぞれについて、対象時刻(t2)における第1変数(xi(t2))を算出し、Xメモリ66に記憶されている直前時刻(t1)における第1変数(xi(t1))を書き換える。これとともに、更新回路69は、N個の要素のそれぞれについて、対象時刻(t2)における第2変数(yi(t2))を算出し、Yメモリ67に記憶されている直前時刻(t1)における第2変数(yi(t1))を書き換える。
制御回路70は、対象時刻(t2)を単位時間(Δt)毎に順次に更新することにより、作用演算回路68および更新回路69に対して単位時間(Δt)毎の第1変数(xi(t))および第2変数(yi(t))を順次に算出させる。
さらに、制御回路70は、1からNまでのインクリメントすることによりインデックス(i)を発生して、作用演算回路68および更新回路69に対して、N個の要素のそれぞれに対応する、対象時刻(t2)における第1変数(xi(t2))および対象時刻(t2)における第2変数(yi(t2))をインデックス順に算出させる。なお、作用演算回路68および更新回路69は、複数個のインデックスに対応する複数個の第1変数(xi(t2))および複数個の第2変数(yi(t2))を並列に算出してもよい。
作用演算回路68は、Jメモリ71と、Hメモリ72と、行列演算回路73と、α関数回路74と、第1加算回路75とを有する。
Jメモリ71は、(N×N)個の結合係数を含むN×Nの行列を記憶する。Ji,jは、行列に含まれるi行j列の結合係数を表す。Ji,jは、組合せ最適化問題を表すイジングモデルにおける、i番目のスピンとj番目のスピンとの結合係数を表す。設定回路64は、演算に先だって、ユーザ等により予め生成された行列をJメモリ71に書き込む。
Hメモリ72は、N個の局所磁場係数を含む配列を記憶する。hiは、配列に含まれるi番目の局所磁場係数を表す。hiは、組合せ最適化問題を表すイジングモデルにおける、i番目のスピンに作用する局所磁場を表す。設定回路64は、演算に先だって、ユーザ等により予め生成された配列をHメモリ72に書き込む。
行列演算回路73は、Xメモリ66から、直前時刻(t1)におけるN個の第1変数(xj(t1))を取得する。行列演算回路73は、N個の要素のそれぞれについて、Jメモリ71から対象の行に含まれるN個の結合係数Ji,jを取得する。そして、行列演算回路73は、N個の要素のそれぞれについて、直前時刻(t1)におけるN個の第1変数(xj(t1))と、対象の行に含まれるN個の結合係数Ji,jとの積和演算を実行する。
α関数回路74は、N個の要素のそれぞれについて、Hメモリ72から対象の局所磁場係数hiを取得する。α関数回路74は、N個の要素のそれぞれについて、{-hiα(t1)}の演算を実行する。α(t)は、予め設定された関数である。
第1加算回路75は、N個の要素のそれぞれについて、行列演算回路73による積和演算結果と、α関数回路74による演算結果とを加算する。この演算により、第1加算回路75は、N個の要素のそれぞれについて、式(43)で表される直前時刻(t
1)における更新値(z
i(t
1))を出力する。
更新回路69は、第1乗算回路79と、P関数回路80と、第2乗算回路81と、第2加算回路82と、第3乗算回路83と、第3加算回路84と、制約前Yメモリ85と、第4乗算回路86と、第4加算回路87と、平均化回路90と、判定回路91と、X制約回路92と、Y制約回路93とを有する。
第1乗算回路79は、N個の要素のそれぞれについて、直前時刻(t1)における更新値(zi(t1))に係数である-cを乗算する。P関数回路80は、N個の要素のそれぞれについて、{-D+p(t1)}の演算を実行する。第2乗算回路81は、N個の要素のそれぞれについて、Xメモリ66から直前時刻(t1)における第1変数(xi(t1))を取得する。そして、第2乗算回路81は、N個の要素のそれぞれについて、直前時刻(t1)における第1変数(xi(t1))と、P関数回路80の演算結果とを乗算する。
第2加算回路82は、N個の要素のそれぞれについて、第1乗算回路79の演算結果と、第2乗算回路81の演算結果とを加算する。第3乗算回路83は、第2加算回路82の演算結果に単位時間であるΔtを乗算する。
第3加算回路84は、N個の要素のそれぞれについて、Yメモリ67から直前時刻(t
1)における第2変数(y
i(t
1))を取得する。第3加算回路84は、N個の要素のそれぞれについて、直前時刻(t
1)における第2変数(y
i(t
1))と、第3乗算回路83の演算結果とを加算する。この演算により、第3加算回路84は、N個の要素のそれぞれについて、式(44)で表される対象時刻(t
2)における第2変数(y
i(t
2))を出力する。
そして、第3加算回路84は、算出したN個の要素のそれぞれの対象時刻(t2)における第2変数(yi(t2))を、制約前Yメモリ85に書き込む。制約前Yメモリ85は、Y制約回路93による制約前の、対象時刻(t2)におけるN個の第2変数(yi(t2))を記憶する。
第4乗算回路86は、N個の要素のそれぞれについて、制約前Yメモリ85から対象時刻(t2)における、制約前の第2変数(yi(t2))を取得する。第4乗算回路86は、N個の要素のそれぞれについて、対象時刻(t2)における制約前の第2変数(yi(t2))に、{DΔt}を乗算する。
第4加算回路87は、N個の要素のそれぞれについて、Xメモリ66から直前時刻(t
1)における第1変数(x
i(t
1))を取得する。第4加算回路87は、N個の要素のそれぞれについて、直前時刻(t
1)における第1変数(x
i(t
1))と、第4乗算回路86の演算結果とを加算する。この演算により、第4加算回路87は、N個の要素のそれぞれについて、式(45)で表される対象時刻(t
2)における第1変数(x
i(t
2))を出力する。
平均化回路90は、Xメモリ66に記憶された直前時刻(t1)におけるN個の第1変数(x1(t1)~xN(t1))の大きさの平均を表す指標値(xave)を算出する。例えば、指標値(xave)は、直前時刻におけるN個の第1変数(x1(t)~xN(t))の二乗平均平方根または絶対値平均である。
判定回路91は、N個の要素のそれぞれについて、第4加算回路87により算出された対象時刻(t2)における第1変数の絶対値(|xi(t2)|)が、予め定められた第2値より大きいか否かを判断する。例えば、第2値は、+1である。判定回路91は、対象時刻(t2)における第1変数の絶対値(|xi(t2)|)が第2値より大きい場合、X制約回路92およびY制約回路93にイネーブル信号(EN)を与える。
X制約回路92は、N個の要素のそれぞれについて、第4加算回路87により算出された対象時刻(t2)における第1変数(xi(t2))を受け取る。X制約回路92は、N個の要素のそれぞれについて、判定回路91からイネーブル信号(EN)を受け取らなかった場合、第4加算回路87により算出された対象時刻(t2)における第1変数(xi(t2))を、そのままXメモリ66に書き込む。
X制約回路92は、N個の要素のそれぞれについて、判定回路91からイネーブル信号(EN)を受け取った場合、第4加算回路87により算出された対象時刻(t2)における第1変数(xi(t2))に対して制約処理を実行し、制約処理をした後の第1変数(xi(t2))をXメモリ66に書き込む。
ここで、X制約回路92は、N個の要素のそれぞれについて、制約処理として、対象時刻(t2)における第1変数(xi(t2))を、符号を変更せずに、絶対値を0以上、第2値(例えば+1)以下の値に変更する。
例えば、N個の要素のそれぞれについて、X制約回路92は、第1変数の絶対値(|xi(t2)|)を、予め定められた値、または、0以上第2値(例えば+1)以下の範囲における所定区間で一様な確率で発生する乱数に応じた値に設定する。この場合、所定区間は、0以上、第2値以下の範囲の内側であれば、どのような範囲であってもよい。
また、N個の要素のそれぞれについて、X制約回路92は、象時刻(t2)における第1変数の絶対値(|xi(t2)|)を、平均化回路90で算出した指標値に変更してもよい。
また、N個の要素のそれぞれについて、X制約回路92は、対象時刻(t2)における第1変数の絶対値(|xi(t2)|)を、平均化回路90で算出した指標値以上、第2値(例えば+1)以下において、乱数により定まる値に変更してもよい。この場合、乱数により定まる値は、例えば、指標値と第2値(例えば+1)以下の範囲内の所定区間で一様な確率で発生する乱数である。
また、N個の要素のそれぞれについて、X制約回路92は、対象時刻(t2)における第1変数の絶対値(|xi(t2)|)を、0から第2値以下の間で初期時刻から終了時刻まで時間経過に従って増加する増加係数に変更してもよい。この場合、増加関数は、例えば、初期時刻で0、終了時刻で第2値となる時刻を変数とする一次関数、または、この一次関数の平方根である。また、N個の要素のそれぞれについて、X制約回路92は、対象時刻(t2)における第1変数の絶対値(|xi(t2)|)を、増加係数以上、第2値(例えば+1)以下における、乱数により定まる値に変更してもよい。この場合、乱数により定まる値は、例えば、増加係数と第2値(例えば+1)との間で一様な確率で発生する乱数である。
Y制約回路93は、N個の要素のそれぞれについて、制約前Yメモリ85から、対象時刻(t2)における第2変数(yi(t2))を取得する。Y制約回路93は、N個の要素のそれぞれについて、判定回路91からイネーブル信号(EN)を受け取らなかった場合、制約前Yメモリ85から取得した対象時刻(t2)における第2変数(yi(t2))を、そのままYメモリ67に書き込む。
Y制約回路93は、N個の要素のそれぞれについて、判定回路91からイネーブル信号(EN)を受け取った場合、第4加算回路87により算出された対象時刻(t2)における第2変数(yi(t2))に対して制約処理を実行し、制約処理をした後の第2変数(yi(t2))をYメモリ67に書き込む。
ここで、Y制約回路93は、N個の要素のそれぞれについて、制約処理として、対象時刻(t2)における第2変数(yi(t2))を、第2変数(yi(t2))に乱数を乗じた値、0、予め定められた値、または、乱数に応じた値に変更する。なお、この場合において、Y制約回路93は、変更後の値が、所定の範囲内の値とするように処理を行う。例えば、Y制約回路93は、変更後の値が、-0.1以上、+0.1以下の範囲内となるように処理を行ってもよい。
以上により、演算回路61は、式(9)の第1アルゴリズムを実行して、終了時刻(T)におけるN個の第1変数(xi(T))およびN個の第2変数(yi(T))を算出することができる。さらに、演算回路61は、更新によって、対象時刻(t2)における第1変数(xi(t1))の絶対値が第2値より大きくなった場合、対象時刻(t2)における第1変数(xi(t1))を、符号を保持したまま、絶対値を0以上第2値以下の値に変更することができる。
なお、上述した演算回路61は、単位時間毎に、対象時刻(t2)における第2変数(yi(t2))を算出した後に、対象時刻(t2)における第1変数(xi(t2))を算出する。これに代えて、演算回路61は、単位時間毎に、対象時刻(t2)における第1変数(xi(t2))を算出した後に、対象時刻(t2)における第2変数(yi(t2))を算出してもよい。
この場合、制約前Yメモリ85は、Y制約回路93による制約前の、直前時刻(t
1)におけるN個の第1変数(x
i(t
1))を記憶する。そして、第4加算回路87は、N個の要素のそれぞれについて、式(46)で表される対象時刻(t
2)における第1変数(x
i(t
2))を出力する。
行列演算回路73は、Xメモリ66から、対象時刻(t
2)におけるN個の第1変数(x
j(t
2))を取得する。そして、第1加算回路75は、N個の要素のそれぞれについて、式(47)で表される対象時刻(t
2)における更新値(z
i(t
2))を出力する。
また、第1乗算回路79は、N個の要素のそれぞれについて、対象時刻(t2)における更新値(zi(t2))に-cを乗算する。P関数回路80は、N個の要素のそれぞれについて、{-D+p(t2)}の演算を実行する。第2乗算回路81は、N個の要素のそれぞれについて、直前時刻(t2)における第1変数(xi(t2))と、P関数回路80の演算結果とを乗算する。
そして、第3加算回路84は、N個の要素のそれぞれについて、式(48)で表される対象時刻(t
2)における第2変数(y
i(t
2))を出力する。
このような処理によっても、演算回路61は、式(9)の第1アルゴリズムを実行して、終了時刻(T)におけるN個の第1変数(xi(T))およびN個の第2変数(yi(T))を算出することができる。さらに、演算回路61は、更新によって、対象時刻(t2)における第1変数(xi(t2))の絶対値が第2値より大きくなった場合、対象時刻(t2)における第1変数(xi(t2))を、符号を保持したまま、絶対値を0以上第2値以下の値に変更することができる。
(第2アルゴリズムの処理を実行する演算回路61)
図26は、式(11)の第2アルゴリズムを実行する演算回路61のブロック構成を示す図である。演算回路61は、式(11)の第2アルゴリズムを実行する場合、図25に示すように構成される。なお、図26に示す演算回路61は、図25に示した構成とほぼ同一の構成を有するので、同一の機能の構成要素については同一の符号を付けて詳細な説明を省略する。
式(11)の第2アルゴリズムを実行する演算回路61は、図25に示した構成と比較して、符号化回路96をさらに有する。
符号化回路96は、Xメモリ66から直前時刻(t1)におけるN個の第1変数(xj(t1))のそれぞれを取得する。符号化回路96は、直前時刻(t1)におけるN個の第1変数(xj(t1))のそれぞれの符号(-1または+1)を抽出し、N個の第1変数の符号(sj(t1))を出力する。
行列演算回路73は、直前時刻(t1)におけるN個の第1変数(xj(t1))に代えて、符号化回路96から、N個の第1変数の符号(sj(t1))を取得する。そして、行列演算回路73は、N個の要素のそれぞれについて、直前時刻(t1)におけるN個の第1変数の符号(sj(t1))と、対象の行に含まれるN個の結合係数Ji,jとの積和演算を実行する。
以上により、演算回路61は、式(11)の第2アルゴリズムを実行して、終了時刻(T)におけるN個の第1変数(xi(T))およびN個の第2変数(yi(T))を算出することができる。
なお、演算回路61が、単位時間毎に、対象時刻(t2)における第1変数(xi(t2))を算出した後に、対象時刻(t2)における第2変数(yi(t2))を算出する場合、符号化回路96は、Xメモリ66から対象時刻(t2)におけるN個の第1変数(xj(t2))のそれぞれを取得する。符号化回路96は、対象時刻(t2)におけるN個の第1変数(xj(t2))のそれぞれの符号(-1または+1)を抽出し、N個の第1変数の符号(sj(t2))を出力する。この場合、行列演算回路73は、対象時刻(t2)におけるN個の第1変数(xj(t2))に代えて、符号化回路96から、N個の第1変数の符号(sj(t2))を取得する。そして、行列演算回路73は、N個の要素のそれぞれについて、対象時刻(t2)におけるN個の第1変数の符号(sj(t2))と、対象の行に含まれるN個の結合係数Ji,jとの積和演算を実行する。このような処理をしても、演算回路61は、式(11)の第2アルゴリズムを実行して、終了時刻(T)におけるN個の第1変数(xi(T))およびN個の第2変数(yi(T))を算出することができる。
(第3アルゴリズムの処理を実行する演算回路61)
図27は、式(13)の第3アルゴリズムを実行する演算回路61のブロック構成を示す図である。演算回路61は、式(13)の第3アルゴリズムを実行する場合、図27に示すように構成される。なお、図27に示す演算回路61は、図25に示した構成とほぼ同一の構成を有するので、同一の機能の構成要素については同一の符号を付けて詳細な説明を省略する。
式(13)の第3アルゴリズムを実行する演算回路61は、図25に示した構成と比較して、更新回路69の内部構成が異なる。式(13)の第3アルゴリズムを実行する更新回路69は、図25に示した構成と比較して、符号化回路96と、G関数回路97とをさらに有する。
符号化回路96は、N個の要素のそれぞれについて、作用演算回路68から、直前時刻(t1)における更新値(zi(t1))を取得する。符号化回路96は、N個の要素のそれぞれについて、直前時刻(t1)における更新値(zi(t1))から符号(-1または+1)を抽出し、直前時刻(t1)における更新値の符号{sgn(zi(t1))}を出力する。
G関数回路97は、N個の要素のそれぞれについて、予め定められた関数g(t)の演算を実行する。具体的には、G関数回路97は、g(t1)={D-p(t1)}√(p(t1))の演算を実行する。
そして、第1乗算回路79は、N個の要素のそれぞれについて、符号化回路96から出力された直前時刻(t1)における更新値の符号{sgn(zi(t1))}と、G関数回路97の演算結果とを乗算する。
以上により、演算回路61は、式(13)の第3アルゴリズムを実行して、終了時刻(T)におけるN個の第1変数(xi(T))およびN個の第2変数(yi(T))を算出することができる。
なお、演算回路61が、単位時間毎に、対象時刻(t2)における第1変数(xi(t2))を算出した後に、対象時刻(t2)における第2変数(yi(t2))を算出する場合、符号化回路96は、作用演算回路68から、対象時刻(t2)における更新値(zi(t2))を取得する。符号化回路96は、N個の要素のそれぞれについて、対象時刻(t2)における更新値(zi(t2))から符号(-1または+1)を抽出し、対象時刻(t2)における更新値の符号{sgn(zi(t2))}を出力する。この場合、G関数回路97は、N個の要素のそれぞれについて、g(t2)={D-p(t2)}√(p(t2))の演算を実行する。そして、第1乗算回路79は、N個の要素のそれぞれについて、符号化回路96から出力された対象時刻(t2)における更新値の符号{sgn(zi(t2))}と、G関数回路97の演算結果とを乗算する。
なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。