JP2022048780A - 計算装置、計算方法およびプログラム - Google Patents

Abstract

【課題】組合せ最適化問題の解を実用的な時間内で計算する。【解決手段】計算装置は、更新部と、出力部とを備える。更新部は、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、単位時間毎の第１変数および第２変数を、単位時間毎に順次に且つ第１変数および第２変数を交互に更新する。単位時間毎の算出処理において、更新部は、複数の要素のそれぞれについて、第１変数を第２変数に基づき更新し、第２変数を複数の要素のそれぞれの第１変数に基づき更新し、第１変数が第１値より小さい場合、第１変数を、第１値以上しきい値以下の値に変更し、第１変数が第２値より大きい場合、第１変数をしきい値以上第２値以下の値に変更する。【選択図】図１８

Description

本発明の実施形態は、計算装置、計算方法およびプログラムに関する。

組合せ最適化問題とは、複数の組合せの中から目的に最も適した組合せを選ぶ問題である。組合せ最適化問題は、数学的には、「目的関数」と呼ばれる、複数の離散変数を有する関数を最大化させる問題、または、当該関数を最小化させる問題に帰着される。組合せ最適化問題は、金融、物流、交通、設計、製造、生命科学など各種の分野において普遍的な問題であるが、組合せ数が問題サイズの指数関数のオーダーで増える、いわゆる「組合せ爆発」のため、必ず最適解を求めることができるとは限らない。また、最適解に近い近
似解を得ることすら難しい場合が多い。

各分野における問題を解決し、社会のイノベーションおよび科学技術の進歩を促進するために、組合せ最適化問題の解を実用的な時間内で計算する技術の開発が求められている。

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

A. Lucas, "Ising formulations of many NP problems", Frontiers in Physics 2, 5 (2014), DOI: 10.3389/fphy.2014.00005 M. W. Johnson et al., "Quantum annealing with manufactured spins", Nature 473, 194-198, 11 May. 2011 T. Inagaki et al.,"A coherent Ising machine for 2000-node optimization problems", Science 354, 603 (2016). H. Goto,"Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network", Sci. Rep. 6, 21686 (2016). M. Yamaoka et al., H. IEEE J. Solid-State Circuits 51, 303 (2016). S. Tsukamoto et al., FUJITSU Sci. Tech. J. 53, 8 (2017). Hayato Goto, Kosuke Tatsumura, Alexander R. Dixon, "Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems",Science Advances, Vol. 5, no. 4, eaav2372, 19 Apr. 2019 EGOR S. TIUNOV, ALEXANDER E. ULANOV, AND A. I. LVOVSKY, "Annealing by simulating the coherent Ising machine",OPTICS EXPRESS 10291, Vol. 27, no. 7, Apr. 2019

本発明の実施形態は、組合せ最適化問題の解を実用的な時間内で計算する計算装置、計算方法およびプログラムを提供する。

実施形態に係る計算装置は、組合せ最適化問題を解く。前記計算装置は、更新部と、出力部とを備える。前記更新部は、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎の前記第１変数および前記第２変数を、前記単位時間毎に順次に且つ前記第１変数および前記第２変数を交互に更新する。前記出力部は、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する。前記複数の要素は、前記組合せ最適化問題を表す複数の離散変数に対応する。前記第１変数および前記第２変数のそれぞれは、実数により表される。前記単位時間毎の算出処理において、前記更新部は、前記複数の要素のそれぞれについて、前記第１変数を前記第２変数に基づき更新し、前記第２変数を前記複数の要素のそれぞれの前記第１変数に基づき更新し、前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値以上、予め定められたしきい値以下の値に変更し、前記第１変数が予め定められた第２値より大きい場合、前記第１変数を、前記しきい値以上、前記第２値以下の値に変更する。前記第２値は、前記第１値より大きい。前記しきい値は、前記第１値より大きく前記第２値より小さい。

情報処理システムの構成例を示した図。 管理サーバの構成例を示したブロック図。 管理サーバの記憶部に保存されるデータの例を示す図。 計算サーバの構成例を示したブロック図。 計算サーバのストレージに保存されるデータの例を示す図。 式（８）のアルゴリズムを使った場合のカット数の平均値を示す図。 式（８）のアルゴリズムを使った場合のカット数の最大値を示す図。 ｘ_ｉの絶対値が１を超うるアルゴリズムにおける分岐の例を示す図。 ｘ_ｉの絶対値が１を超えないアルゴリズムにおける分岐の例を示す図。 第１アルゴリズムを使った場合のカット数の平均値を示す図。 第１アルゴリズムを使った場合のカット数の最大値を示す図。 第２アルゴリズムを使った場合のカット数の平均値を示す図。 第２アルゴリズムを使った場合のカット数の最大値を示す図。 第３アルゴリズムを使った場合のカット数の平均値を示す図。 第３アルゴリズムを使った場合のカット数の最大値を示す図。 ｃ１を適用した第１アルゴリズムを使った場合のカット数の平均値を示す図。 ｃ１を適用した第１アルゴリズムを使った場合のカット数の最大値を示す図。 ｃ１を適用した第２アルゴリズムを使った場合のカット数の平均値を示す図。 ｃ１を適用した第２アルゴリズムを使った場合のカット数の最大値を示す図。 マルチプロセッサ構成の例を概略的に示した図。 ＰＣクラスタで問題を解いた場合のエネルギーの平均値を示す図。 ＰＣクラスタで問題を解いた場合の計算時間の平均値を示す図。 ＧＰＵを使った構成の例を概略的に示した図。 ＧＰＵで問題を解いた場合のエネルギーの平均値を示す図。 ＧＰＵで問題を解いた場合の計算時間の平均値を示す図。 第１アルゴリズムを実行する場合における処理の流れの第１例を示す図。 第１アルゴリズムを実行する場合における処理の流れの第２例を示す図。 第２アルゴリズムを実行する場合における処理の流れの第１例を示す図。 第２アルゴリズムを実行する場合における処理の流れの第２例を示す図。 第３アルゴリズムを実行する場合における処理の流れの第１例を示す図。 第３アルゴリズムを実行する場合における処理の流れの第２例を示す図。 演算装置の構成を示す図。 第１アルゴリズムを実行する演算回路のブロック構成図。 第２アルゴリズムを実行する演算回路のブロック構成図。 第３アルゴリズムを実行する演算回路のブロック構成図。

以下、図面を参照しながら、本発明の実施形態について説明する。また、図面において同一の構成要素は、同じ番号を付し、説明は、適宜省略する。

（システム構成）
図１は、情報処理システム１００の構成例を示したブロック図である。図１の情報処理システム１００は、管理サーバ１と、ネットワーク２と、複数の計算サーバ（情報処理装置）３（３ａ～３ｃ）と、複数のケーブル４（４ａ～４ｃ）と、スイッチ５を備えている。また、図１には、情報処理システム１００と通信可能な情報端末６が示されている。管理サーバ１、複数の計算サーバ３（３ａ～３ｃ）、情報端末６は、ネットワーク２を介して互いにデータ通信をすることができる。ネットワーク２は、例えば、複数のコンピュータネットワークが相互に接続されたインターネットである。ネットワーク２は、通信媒体として有線、無線、または、これらの組み合わせを用いることができる。また、ネットワーク２で使われる通信プロトコルの例としては、ＴＣＰ／ＩＰがあるが、通信プロトコルの種類については特に問わない。

また、複数の計算サーバ３（３ａ～３ｃ）は、それぞれケーブル４（４ａ～４ｃ）を介してスイッチ５に接続されている。複数のケーブル４（４ａ～４ｃ）およびスイッチ５は、計算サーバ間のインターコネクトを形成している。複数の計算サーバ３（３ａ～３ｃ）のそれぞれは、当該インターコネクトを介して互いにデータ通信をすることも可能である。スイッチ５は、例えば、Ｉｎｆｉｎｉｂａｎｄのスイッチであり、ケーブル４ａ～４ｃは、例えば、Ｉｎｆｉｎｉｂａｎｄのケーブルである。ただし、Ｉｎｆｉｎｉｂａｎｄのスイッチ／ケーブルの代わりに、有線ＬＡＮのスイッチ／ケーブルを使ってもよい。ケーブル４ａ～４ｃおよびスイッチ５で使われる通信規格および通信プロトコルについては、特に問わない。情報端末６の例としては、ノートＰＣ、デスクトップＰＣ、スマートフォン、タブレット、車載端末などが挙げられる。

組合せ最適化問題の求解では、並列的な処理および／または処理の分散化を行うことができる。したがって、複数の計算サーバ３（３ａ～３ｃ）のそれぞれおよび／または計算サーバ３（３ａ～３ｃ）のプロセッサは、一部の計算処理の一部のステップを分担して実行してもよいし、異なる変数について同様の計算処理を並列的に実行してもよい。管理サーバ１は、例えば、ユーザによって入力された組合せ最適化問題を各計算サーバ３に処理可能な形式に変換し、計算サーバ３を制御する。そして、管理サーバ１は、各計算サーバ３から計算結果を取得し、集約した計算結果を組合せ最適化問題の解に変換する。これにより、ユーザは、組合せ最適化問題の解を得ることができる。組合せ最適化問題の解は、最適解と、最適解に近い近似解とを含むものとする。

図１には、３台の計算サーバ３（３ａ～３ｃ）が示されているが、これは、情報処理システム１００に含まれる計算サーバ３の台数を限定することを意図していない。また、組合せ最適化問題の求解に使われる計算サーバ３の台数についても特に問わない。例えば、情報処理システム１００に含まれる計算サーバ３は１台であってもよい。また、情報処理システム１００に含まれる複数の計算サーバ３のうち、１台の計算サーバ３を使って組合せ最適化問題の求解を行ってもよい。また、情報処理システム１００に、数百台以上の計算サーバ３が含まれていてもよい。計算サーバ３は、データセンターに設置されたサーバであってもよいし、オフィスに設置されたデスクトップＰＣであってもよい。また、計算サーバ３は異なるローケーションに設置された複数の種類のコンピュータであってもよい。計算サーバ３として使われる情報処理装置の種類については特に問わない。例えば、計算サーバ３は、汎用的なコンピュータであってもよいし、専用の電子回路または、これらの組合せであってもよい。

図２は、管理サーバ１の構成例を示したブロック図である。図２の管理サーバ１は、例えば、中央演算処理装置（ＣＰＵ）とメモリとを含むコンピュータである。管理サーバ１は、プロセッサ１０と、記憶部１４と、通信回路１５と、入力回路１６と、出力回路１７とを備えている。プロセッサ１０、記憶部１４、通信回路１５、入力回路１６、出力回路１７は、互いにバス２０を介して接続されているものとする。プロセッサ１０は、内部の機能構成として、管理部１１と、変換部１２と、制御部１３を含んでいる。

プロセッサ１０は、演算を実行し、管理サーバ１の制御を行う電子回路である。プロセッサ１０として、例えば、ＣＰＵ、マイクロプロセッサ、ＡＳＩＣ、ＦＰＧＡ、ＰＬＤまたはこれらの組み合わせを用いることができる。管理部１１は、ユーザの情報端末６を介して管理サーバ１の操作を行うためのインタフェースを提供する。管理部１１が提供するインタフェースの例としては、ＡＰＩ、ＣＬＩ、ウェブページなどが挙げられる。例えば、ユーザは、管理部１１を介して組合せ最適化問題の情報の入力を行ったり、計算された組合せ最適化問題の解の閲覧および／またはダウンロードを行ったりすることができる。変換部１２は、組合せ最適化問題に関するパラメータを入力して、入力したパラメータを各計算サーバ３が処理可能な形式に変換する。制御部１３は、各計算サーバ３に制御指令を送信する。制御部１３が各計算サーバ３から計算結果を取得した後、変換部１２は、複数の計算結果を集約し、組合せ最適化問題の解に変換し、組合せ最適化問題の解を出力する。

記憶部１４は、管理サーバ１のプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、ＯＳとアプリケーションの両方を含むものとする。記憶部１４は、揮発性メモリ、不揮発性メモリ、またはこれらの組み合わせであってもよい。揮発性メモリの例としては、ＤＲＡＭ、ＳＲＡＭなどがある。不揮発性メモリの例としては、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭ、ＭＲＡＭが挙げられる。また、記憶部１４として、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置を使ってもよい。

通信回路１５は、ネットワーク２に接続された各装置との間でデータの送受信を行う。通信回路１５は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ Ｉｎｔｅｒｆａｃｅ Ｃａｒｄ）である。ただし、通信回路１５は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。入力回路１６は、管理サーバ１へのデータ入力を実現する。入力回路１６は、外部ポートとして、例えば、ＵＳＢ、ＰＣＩ－Ｅｘｐｒｅｓｓなどを備えているものとする。図２の例では、操作装置１８が入力回路１６に接続されている。操作装置１８は、管理サーバ１に情報を入力するための装置である。操作装置１８は、例えば、キーボード、マウス、タッチパネル、音声認識装置などであるが、これに限られない。出力回路１７は、管理サーバ１からのデータ出力を実現する。出力回路１７は、外部ポートとしてＨＤＭＩ（登録商標）、ＤｉｓｐｌａｙＰｏｒｔなどを備えているものとする。図２の例では、表示装置１９が出力回路１７に接続されている。表示装置１９の例としては、ＬＣＤ（液晶ディスプレイ）、有機ＥＬ（有機エレクトロルミネッセンス）ディスプレイ、プロジェクタがあるが、これに限られない。

管理サーバ１の管理者は、操作装置１８および表示装置１９を使って、管理サーバ１のメンテナンスを行うことができる。なお、操作装置１８および表示装置１９は、管理サーバ１に組み込まれたものであってもよい。また、必ず管理サーバ１に操作装置１８および表示装置１９が接続されていなくてもよい。例えば、管理者は、ネットワーク２と通信可能な情報端末を用いて管理サーバ１のメンテナンスを行ってもよい。

図３は、管理サーバ１の記憶部１４に保存されるデータの例を示している。図３の記憶部１４には、問題データ１４Ａと、計算データ１４Ｂと、管理プログラム１４Ｃと、変換プログラム１４Ｄと、制御プログラム１４Ｅとが保存されている。例えば、問題データ１４Ａは、組合せ最適化問題のデータを含む。例えば、計算データ１４Ｂは、各計算サーバ３から収集された計算結果を含む。例えば、管理プログラム１４Ｃは、上述の管理部１１の機能を実現するプログラムである。例えば、変換プログラム１４Ｄは、上述の変換部１２の機能を実現するプログラムである。例えば、制御プログラム１４Ｅは、上述の制御部１３の機能を実現するプログラムである。

図４は、計算サーバ３ａの構成例を示したブロックである。図４には、例示的に計算サーバ３ａの構成が示されている。他の計算サーバ３は、計算サーバ３ａと同様の構成であってもよいし、計算サーバ３ａと異なる構成であってもよい。計算サーバ３ａは、例えば、第１ベクトルと、第２ベクトルと、第３ベクトルの計算を単独で、または、他の計算サーバ３と分担して実行する情報処理装置である。計算サーバ３ａは、さらに第１ベクトルのそれぞれの要素を符号関数で変換した第４ベクトルの計算を行ってもよい。第３ベクトルのそれぞれの要素の値は、例えば、イジングモデルのエネルギー式から導出された式によって求められる。例えば、イジングモデルのエネルギー式を、すべての項に含まれる変数について偏微分した形式の式（基本式とよぶ）に基づいて第３ベクトルのそれぞれの要素を計算することができる。

ここで、第１ベクトルは、変数ｘ_ｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルである。また、第２ベクトルは、変数ｙ_ｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルである。第３ベクトルは、変数ｚ_ｉ（ｉ＝１、２、・・・、Ｎ）を要素とするベクトルである。第４ベクトルは、第１ベクトルの要素を第１値もしくは第１値より大きい第２値のいずれかの値をとる第１関数で変換したベクトルである。上述の符号関数は、第１関数の一例である。なお、変数ｘ_ｉ、ｙ_ｉ、ｚ_ｉの詳細については、後述する。

計算サーバ３ａは、例えば、通信回路３１と、共有メモリ３２と、プロセッサ３３Ａ～３３Ｄと、ストレージ３４と、ホストバスアダプタ３５とを備えている。通信回路３１、共有メモリ３２、プロセッサ３３Ａ～３３Ｄ、ストレージ３４、ホストバスアダプタ３５は、バス３６を介して互いに接続されているものとする。

通信回路３１は、ネットワーク２に接続された各装置との間でデータの送受信を行う。通信回路３１は、例えば、有線ＬＡＮのＮＩＣ（Ｎｅｔｗｏｒｋ Ｉｎｔｅｒｆａｃｅ Ｃａｒｄ）である。ただし、通信回路３１は、無線ＬＡＮなど、その他の種類の通信回路であってもよい。共有メモリ３２は、プロセッサ３３Ａ～３３Ｄからアクセス可能なメモリである。共有メモリ３２の例としては、ＤＲＡＭ、ＳＲＡＭなどの揮発性メモリが挙げられる。ただし、共有メモリ３２として、不揮発性メモリなどその他の種類のメモリが使われてもよい。プロセッサ３３Ａ～３３Ｄは、共有メモリ３２を介してデータの共有を行うことができる。なお、必ず計算サーバ３ａのすべてのメモリが共有メモリとして構成されていなくてもよい。例えば、計算サーバ３ａの一部のメモリは、いずれかのプロセッサのみからアクセスできるローカルメモリとして構成されていてもよい。

プロセッサ３３Ａ～３３Ｄは、計算処理を実行する電子回路である。プロセッサは、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）のいずれであってもよいし、これらの組合せであってもよい。また、プロセッサは、ＣＰＵコアまたはＣＰＵスレッドであってもよい。プロセッサがＣＰＵである場合、計算サーバ３ａが備えるソケット数については、特に問わない。また、プロセッサは、ＰＣＩ ｅｘｐｒｅｓｓなどのバスを介して計算サーバ３ａのその他の構成要素に接続されていてもよい。

図４の例では、計算サーバ３ａが４つのプロセッサを備えている。ただし、１台の計算サーバ３ａが備えているプロセッサの数はこれとは異なっていてもよい。例えば、計算サーバ３ａに実装されているプロセッサの数および／または種類が異なっていてもよい。

作用演算部５１は、解きたい組合せ最適化問題の目的関数を、すべての項に含まれる変数について偏微分した形式の基本式に基づき、第３ベクトルのそれぞれの要素を更新するように構成されている。ここで、基本式の変数は、第１ベクトルの要素、または、第１ベクトルの要素を第１値もしくは第１値より大きい第２値のいずれかの値をとる第１関数で変換した第４ベクトルの要素である。更新部５０は、例えば、第１ベクトルの要素に、第２ベクトルの対応する要素、または、第２ベクトルの対応する要素に重み付けした値を加算することによって第１ベクトルの要素を更新し、値が第１値より小さい第１ベクトルの要素を第１値以上しきい値以下の何れかの値に変更し、値が第２値より大きい第１ベクトルの要素をしきい値以上第２値以下の値に変更し、第２ベクトルの要素に、更新回数に応じて単調増加または単調減少する第１係数と第１ベクトルの対応する要素との積を重み付けした値、ならびに、第３ベクトルの対応する要素を重み付けした値を加算することによって第２ベクトルの要素を更新するように構成されている。なお、しきい値は、第１値と第２値との間の値である。目的関数として、例えば、イジングモデルのエネルギー式を使うことができる。ここで、イジングモデルは、多体相互作用を有するものであってもよい。また、第１値として－１、第２値として＋１、しきい値として０を使うことができる。ただし、しきい値、第１値および／または第２値はその他の値であってもよい。

図４の例では、プロセッサ３３Ａ～３３Ｃが更新部５０に相当しており、プロセッサ３３Ｄが作用演算部５１に相当している。ただし、図４に示した更新部５０／作用演算部５１とプロセッサとの対応関係は、一例にしかすぎない。したがって、更新部５０／作用演算部５１とプロセッサとの対応関係はこれとは異なっていてもよい。また、更新部５０／作用演算部５１に割り当てられるプロセッサ数については、特に限定しない。後述するように、同一のプロセッサが更新部５０および作用演算部５１の役割を兼ね備えていてもよい。計算サーバ３ａに複数の種類のプロセッサ（例えば、ＣＰＵ、ＧＰＵ、ＦＰＧＡ）が実装されている場合には、異なる種類のプロセッサを更新部５０および作用演算部５１に割り当ててもよい。

ストレージ３４は、計算サーバ３ａのプログラム、プログラムの実行に必要なデータ、プログラムによって生成されたデータを含む各種のデータを記憶する。ここで、プログラムは、ＯＳとアプリケーションの両方を含むものとする。ストレージ３４は、揮発性メモリ、不揮発性メモリ、またはこれらの組み合わせであってもよい。揮発性メモリの例としては、ＤＲＡＭ、ＳＲＡＭなどがある。不揮発性メモリの例としては、ＮＡＮＤフラッシュメモリ、ＮＯＲフラッシュメモリ、ＲｅＲＡＭ、ＭＲＡＭが挙げられる。また、ストレージ３４として、ハードディスク、光ディスク、磁気テープまたは外部の記憶装置が使われてもよい。

ホストバスアダプタ３５は、計算サーバ３間のデータ通信を実現する。ホストバスアダプタ３５は、ケーブル４ａを介してスイッチ５に接続されている。ホストバスアダプタ３５は、例えば、ＨＣＡ（Ｈｏｓｔ Ｃｈａｎｎｅｌ Ａｄａｐｔｏｒ）である。ホストバスアダプタ３５、ケーブル４ａ、スイッチ５で高スループットを実現可能なインターコネクトを形成することにより、並列的な計算処理の速度を向上させることができる。

図５は、計算サーバ３のストレージに保存されるデータの例を示している。図５のストレージ３４には、計算データ３４Ａと、計算プログラム３４Ｂと、制御プログラム３４Ｃとが保存されている。計算データ３４Ａは、計算サーバ３ａの計算途中のデータまたは計算結果を含んでいる。なお、計算データ３４Ａの少なくとも一部は、共有メモリ３２、プロセッサのキャッシュ、プロセッサのレジスタなど、異なる記憶階層に保存されていてもよい。計算プログラム３４Ｂは、所定のアルゴリズムに基づき、各プロセッサにおける計算処理および、共有メモリ３２およびストレージ３４へのデータの保存処理を実現するプログラムである。制御プログラム３４Ｃは、管理サーバ１の制御部１３から送信された指令に基づき、計算サーバ３ａを制御し、計算サーバ３ａの計算結果を管理サーバ１に送信するプログラムである。

（組合せ最適化問題）
次に、組合せ最適化問題の求解に関連する技術について説明する。組合せ最適化問題を解くために使われる情報処理装置の一例として、イジングマシンが挙げられる。イジングマシンとは、イジングモデルの基底状態のエネルギーを計算する情報処理装置のことをいう。これまで、イジングモデルは、主に強磁性体や相転移現象のモデルとして使われることが多かった。しかし、近年、イジングモデルは、組み合わせ最適化問題を解くためのモデルとしての利用が増えている。下記の式（１）は、イジングモデルのエネルギーを示している。

ここで、ｓ_ｉ、ｓ_ｊはスピンである、スピンは、＋１または－１のいずれかの値をとる２値変数である。Ｎは、スピンの数である。ｈ_ｉは、各スピンに作用する局所磁場である。Ｊは、スピン間における結合係数の行列である。行列Ｊは、対角成分が０である実対称行列となっている。したがって、Ｊ_ｉｊは行列Ｊのｉ行ｊ列の要素を示している。なお、式（１）のイジングモデルは、スピンについての２次式となっているが、スピンの３次以上の項を含む拡張されたイジングモデル（多体相互作用を有するイジングモデル）を使ってもよい。多体相互作用を有するイジングモデルについては、後述する。

式（１）のイジングモデルを使うと、エネルギーＥ_{Ｉｓｉｎｇ}を目的関数とし、エネルギーＥ_{Ｉｓｉｎｇ}を可能な限り小さくする解を計算することができる。イジングモデルの解は、スピンのベクトル（ｓ_１、ｓ_２、・・・、ｓ_Ｎ）の形式で表される。特に、エネルギーＥ_{Ｉｓｉｎｇ}が最小値となるベクトル（ｓ_１、ｓ_２、・・・、ｓ_Ｎ）は、最適解とよばれる。ただし、計算されるイジングモデルの解は、必ず厳密な最適解でなくてもよい。以降では、イジングモデルを使ってエネルギーＥ_{Ｉｓｉｎｇ}が可能な限り小さくなる近似解（すなわち、目的関数の値が可能な限り最適値に近くなる近似解）を求める問題をイジング問題とよぶものとする。

式（１）のｓ_ｉは、スピンを表す２値変数であるため、式（１＋ｓ_ｉ）／２を使うことにより、組合せ最適化問題で使われる離散変数（ビット）との変換を容易に行うことができる。したがって、組合せ最適化問題をイジング問題に変換し、イジングマシンに計算を行わせることにより、組合せ最適化問題の解を求めることが可能である。０または１のいずれかの値をとる離散変数（ビット）を変数とする２次の目的関数を最小化する解を求める問題は、ＱＵＢＯ（Ｑｕａｄｒａｔｉｃ Ｕｎｃｏｎｓｔｒａｉｎｅｄ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ、制約なし２値変数２次最適化）問題とよばれる。式（１）で表されるイジング問題は、ＱＵＢＯ問題と等価であるといえる。

例えば、量子アニーラ、コヒーレントイジングマシン、量子分岐マシンなどがイジングマシンのハードウェア実装として提案されている。量子アニーラは、超伝導回路を使って量子アニーリングを実現する。コヒーレントイジングマシンは、光パラメトリック発振器で形成されたネットワークの発振現象を利用する。量子分岐マシンは、カー効果を有するパラメトリック発振器のネットワークにおける量子力学的な分岐現象を利用する。これらのハードウェア実装は、計算時間の大幅な短縮を実現する可能性がある一方、大規模化や安定的な運用が難しいという課題もある。

そこで、広く普及しているデジタルコンピュータを使ってイジング問題の求解を行うことも可能である。デジタルコンピュータは、上述の物理的現象を使ったハードウェア実装と比べ、大規模化と安定運用が容易である。デジタルコンピュータでイジング問題の求解を行うためのアルゴリズムの一例として、シミュレーティッドアニーリング（ＳＡ）が挙げられる。シミュレーティッドアニーリングをより高速に実行する技術の開発が行われている。ただし、一般のシミュレーティッドアニーリングはそれぞれの変数が逐次更新される逐次更新アルゴリズムであるため、並列化による計算処理の高速化は難しい。

（シミュレーテッド分岐アルゴリズム）
上述の課題を踏まえ、デジタルコンピュータにおける並列的な計算によって、規模の大きい組合せ最適化問題の求解を高速に行うことが可能なシミュレーテッド分岐アルゴリズムが提案されている（例えば、上述の非特許文献７）。以降では、シミュレーテッド分岐アルゴリズムを使って組合せ最適化問題を解く情報処理装置および電子回路について説明する。

はじめに、シミュレーテッド分岐アルゴリズムの概要について述べる。シミュレーテッド分岐アルゴリズムでは、それぞれがＮ個の要素に対応する２つの変数ｘ_ｉおよび変数ｙ_ｉを用いる。なお、変数ｘ_ｉを第１変数、変数ｙ_ｉを第２変数と呼ぶ場合もある。ここで、シミュレーテッド分岐アルゴリズムにおいて、Ｎ個の要素のそれぞれは、仮想的な粒子を表している。Ｎ個の要素は、最適化問題を表すイジングモデルにおけるＮ個のスピンに対応する。第１変数ｘ_ｉは、Ｎ個の粒子のうちのｉ番目の粒子の位置を表す。第２変数ｙ_ｉは、ｉ番目の粒子の運動量を表す。Ｎは、イジングモデルに含まれるスピンの数を表し、２以上の整数である。ｉは、１以上、Ｎ以下の任意の整数を表し、スピンを特定するインデックスを表す。そして、シミュレーテッド分岐アルゴリズムでは、それぞれＮ個ある２つの変数ｘ_ｉ，ｙ_ｉ（ｉ＝１、２、・・・、Ｎ）について、下記の式（２）の連立常微分方程式を数値的に解く。変数ｘ_ｉ，ｙ_ｉは、いずれも、実数で表される連続変数である。

ここで、Ｈは、下記の式（３）のハミルトニアンである。係数Ｄは、予め定められた定数であり、離調（ｄｅｔｕｎｉｎｇ）に相当する。係数ｐ（ｔ）は、ポンピング振幅（ｐｕｍｐｉｎｇ ａｍｐｌｉｔｕｄｅ）に相当し、シミュレーテッド分岐アルゴリズムの計算時に更新回数に応じて値が単調増加する。ｔは、時刻を表す変数である。係数ｐ（ｔ）の初期値は０に設定されていてもよい。係数ｐ（ｔ）は、第１係数に相当する。係数Ｋは、正のカー係数（Ｋｅｒｒ ｃｏｅｆｆｉｃｉｅｎｔ）に相当する。ｆ_ｉは、下記の式（４）で表される外力である。式（４）のｚ_ｉは、式（３）の中のエネルギーＥ_{Ｉｓｉｎｇ}に対応する項の括弧内を変数ｘ_ｉで偏微分した式となっている。

ここで、係数ｃとして、定数係数を使うことができる。この場合、係数ｃの値を、シミュレーテッド分岐アルゴリズムによる計算を実行する前に決める必要がある。例えば、計算の精度を得るために、係数ｃをＪ^（２）行列の最大固有値の逆数に近い値に設定することができる。例えば、ｃ＝０．５Ｄ√（Ｎ／２ｎ）という値を使うことができる。ここで、ｎは、組合せ最適化問題に係るグラフのエッジ数である。また、α（ｔ）は、ｐ（ｔ）とともに増加する係数である。例えば、α（ｔ）として、√（ｐ（ｔ））を使うことができる。

なお、シミュレーテッド分岐アルゴリズムを使うことにより、３次以上の目的関数を有する組合せ最適化問題を解くことも可能である。２値変数を変数とする３次以上の目的関数を最小化する変数の組合せを求める問題は、ＨＯＢＯ（Ｈｉｇｈｅｒ Ｏｒｄｅｒ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）問題とよばれる。ＨＯＢＯ問題を扱う場合、高次へ拡張されたイジングモデルにおけるエネルギー式として、下記の式（５）を使うことができる。

ここで、Ｊ^（ｎ）はｎ階テンソルであり、式（１）の局所磁場ｈ_ｉと結合係数の行列Ｊを一般化させたものである。例えば、テンソルＪ^（１）は、局所磁場ｈ_ｉのベクトル（第６ベクトルとよぶ）に相当する。ｎ階テンソルＪ^（ｎ）では、複数の添え字に同じ値があるとき、要素の値は０となる。式（５）では、３次の項までが示されているが、それより高次の項も式（５）と同様に定義することができる。式（５）は多体相互作用を含むイジングモデルのエネルギーに相当している。

なお、ＱＵＢＯと、ＨＯＢＯはいずれも、制約なし多項式２値変数最適化（ＰＵＢＯ：Ｐｏｌｙｎｏｍｉａｌ Ｕｎｃｏｎｓｔｒａｉｎｅｄ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）の１種であるといえる。すなわち、ＰＵＢＯのうち、２次の目的関数を有する組合せ最適化問題は、ＱＵＢＯである。また、ＰＵＢＯのうち、３次以上の目的関数を有する組合せ最適化問題は、ＨＯＢＯであるといえる。

シミュレーテッド分岐アルゴリズムを使ってＨＯＢＯ問題を解く場合、上述の式（３）のハミルトニアンＨを下記の式（６）に、上述の式（４）の外力ｆ_ｉを下記の式（７）にそれぞれ置き換えればよい。

例えば、（７）の２番目の式ｚ_ｉを用いて、第３ベクトルのそれぞれの要素を計算することができる。この式は、（６）の２番目の式を、すべての項に含まれる変数ｘ_ｉについて偏微分した形式をとっている。また、第１ベクトルの要素を変数としている。このように、ハミルトニアンが多体相互作用（３階以上のテンソル）の項を含んでもよい。また、ハミルトニアンとして、多体相互作用（３階以上のテンソル）の項を含まないものを使ってもよい。（７）の２番目の式ｚ_ｉは、ハミルトニアンの中のイジングモデルのエネルギーに対応する項から導出された基本式の一例である。すなわち、第１値は－１、第２値は１であってもよく、目的関数は、イジングモデルのエネルギー式に相当する項を含んでいてもよい。この場合、目的関数は、多体相互作用の項を含んでいてもよい。

シミュレーテッド分岐アルゴリズムでは、ｐ（ｔ）の値を初期値（例えば、０）から所定の値まで増加させた後における変数ｘ_ｉの符号に基づき、スピンｓ_ｉの値を求めることができる。例えば、ｘ_ｉ＞０のときｓｇｎ（ｘ_ｉ）＝１、ｘ_ｉ＜０のときｓｇｎ（ｘ_ｉ）＝－１となる符号関数を使うと、ｐ（ｔ）の値が所定の値まで増加したとき、変数ｘ_ｉを符号関数で変換することによってスピンｓ_ｉの値を求めることができる。符号関数として、例えば、ｘ_ｉ≠０のときに、ｓｇｎ（ｘ_ｉ）＝ｘ_ｉ／｜ｘ_ｉ｜、ｘ_ｉ＝０のときにｓｇｎ（ｘ_ｉ）＝１または－１になる関数を使うことができる。すなわち、更新部５０は、値が第１値と第２値の間にある第３値より小さい第１ベクトルの要素を第１値に変換し、値が第３値より大きい第１ベクトルの要素を第２値に変換することによって、組合せ最適化問題の解を求めるように構成されていてもよい。例えば、更新部５０は、正値である第１ベクトルの要素を＋１に変換し、負値である第１ベクトルを－１に変換することによって、組合せ最適化問題の解を求めるように構成されていてもよい。更新部５０が組合せ最適化問題の解（例えば、イジングモデルのスピンｓ_ｉ）を求めるタイミングについては、特に問わない。例えば、更新部５０は、第１ベクトル、第２ベクトル、第３ベクトルの更新回数または第１係数ｐの値がしきい値より大きくなったときに組合せ最適化問題の解を求めるように構成されていてもよい。イジング問題を解く場合、組合せ最適化問題の解は、イジングモデルのスピンｓ_ｉに相当する。

（シミュレーテッド分岐アルゴリズムの演算）
例えば、シンプレクティック・オイラー法を使うと、式（２）、（３）、（４）または、式（２）、（６）、（７）によって与えられる微分方程式を解くことができる。下記の式（８）に示されているように、シンプレクティック・オイラー法を使う場合、微分方程式が離散的な漸化式に書き換えられる。

ここで、ｔは、時刻であり、Δｔは、時間ステップ（単位時間、時間刻み幅）である。式（８）の非線形項Ｋｘ^２ _ｉ（ｔ＋Δｔ）は、計算中に変数ｘ_ｉが発散するのを防止する。

計算サーバ３では、式（８）のアルゴリズムに基づき、それぞれＮ個ある２つの変数ｘ_ｉ，ｙ_ｉ（ｉ＝１、２、・・・、Ｎ）を更新してもよい。すなわち、計算サーバ３が更新するデータには、変数ｘ_ｉ（ｉ＝１、２、・・・、Ｎ）を要素とする第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）と、変数ｙ_ｉ（ｉ＝１、２、・・・、Ｎ）を要素とする第２ベクトル（ｙ_１，ｙ_２，・・・，ｙ_Ｎ）と、変数ｚ_ｉ（ｉ＝１、２、・・・、Ｎ）を要素とする第３ベクトル（ｚ_１，ｚ_２，・・・，ｚ_Ｎ）が含まれていてもよい。計算サーバ３は、式（８）のアルゴリズムに基づき、第３ベクトルのそれぞれの要素ｚ_ｉ（ｉ＝１、２、・・・、Ｎ）と、第１のベクトルのそれぞれの要素ｘ_ｉ（ｉ＝１、２、・・・、Ｎ）と、第２ベクトルのそれぞれの要素ｙ_ｉ（ｉ＝１、２、・・・、Ｎ）を更新することができる。

式（８）を参照すると、外力項ｆ_ｉに含まれる、行列またはテンソルの積和演算以外は、１種類の添え字（ｉ）しか現れていないことがわかる。このため、式（８）のうち、１種類の添え字（ｉ）しか現れていない部分の演算を並列化することによって、計算時間を短縮することができる。

なお、式（８）では、微分方程式との対応関係を示すために、時刻ｔおよび時間ステップΔｔが使われている。ただし、実際にシンプレクティック・オイラー法をソフトウェアまたはハードウェアに実装する際は、必ず明示的なパラメータとして時刻ｔおよび時間ステップΔｔが含まれていなくてもよい。例えば、時間ステップΔｔを１とすれば、実装時のアルゴリズムから時間ステップΔｔを除去することが可能である。アルゴリズムを実装する際に、明示的なパラメータとして時刻ｔを含めない場合には、式（８）において、ｘ_ｉ（ｔ＋Δｔ）をｘ_ｉ（ｔ）の更新後の値であると解釈すればよい。すなわち、上述の式（８）および以降の各式における“ｔ”は、更新前の変数の値、“ｔ＋Δｔ”は、更新後の変数の値を示すものとする。

次に、シミュレーテッド分岐アルゴリズムをシンプレクティック・オイラー法によってデジタルコンピュータに実装し、組合せ最適化問題を解いたときの結果を示す。以降では、最大カット問題のベンチマークセット（Ｇ－ｓｅｔ）のＧ２２を１０００回解いた場合におけるカット数の平均値と最大値を示す。最大カット問題とは、分割時にカットされるエッジの重みの合計値を最大化するよう、重み付きグラフのノードを２つのグループに分割する問題である。最大カット問題は、組合せ最適化問題の一種である。

図６Ａおよび図６Ｂは、上述の式（８）のアルゴリズムを使ったときにおける結果を示している。時間ステップをΔｔ＝０．５、合計時間ステップ数を１００、１０００、１００００、１０００００として計算を行っている。係数について、Ｄ＝Ｋ＝１、ｃ=０．５Ｄ√（Ｎ／２ｎ）を使った。なお、ｎにはＧ２２のグラフのエッジ数、１９９９０が代入される。時間ステップ数の増加に応じて係数ｐ（第１係数）の値を０から１に線形に増加させた。また、変数ｘ_ｉの初期値に０を設定し、変数ｙ_ｉの初期値として［－０．１，０．１］の範囲の擬似乱数を設定した。

図６Ａは、カット数の平均値を示している。一方、図６Ｂは、カット数の最大値を示している。図６Ａおよび図６Ｂのいずれのグラフにおいても、縦軸はカット数に、横軸は時間ステップ数にそれぞれ対応している。図６Ａおよび図６Ｂの両グラフにある水平方向の破線Ｃｍａｘは、Ｇ２２で知られている最大カット数１３３５９を示している。カット数が破線Ｃｍａｘに近いほど、最適解に近い結果が得られているといえる。図６Ａおよび図６Ｂを参照すると、合計時間ステップ数が大きくなっても最大カット数に到達していない。

（アルゴリズムの改良）
図７は、式（８）のアルゴリズムの分岐現象を表す図である。シミュレーテッド分岐アルゴリズムをシンプレクティック・オイラー法によって解いた場合、系のパラメータの変化に伴い、安定運動状態が１個から複数個へと分岐する。式（８）のアルゴリズムの分岐現象では、変数ｘ_ｉが１より大きい領域、または、－１より小さい領域まで広がっている。

図８は、改良したアルゴリズムの分岐現象を表す図である。そこで、以下のように、式（８）のアルゴリズムを改良した。具体的には、図８に示すように、更新によって変数ｘ_ｉの絶対値が１より大きくなった場合、変数ｘ_ｉは、符号が変更されずに、絶対値を０以上１以下の値に変更される。例えば、更新によってｘ_ｉ＞１となった場合、変数ｘ_ｉの値はｗに設定される。ここで、ｗは、０以上１以下の値である。また、更新によってｘ_ｉ＜－１となった場合、変数ｘ_ｉの値は－ｗに設定される。これにより、変数ｘ_ｉの符号を保ちつつ、変数ｘ_ｉを常に－１以上、１以下の範囲に保つことができる。

例えば、ｗは、０以上１以下の予め定められた値であってもよい。また、ｗは、０以上１以下の範囲内の所定区間［ｗ_１，ｗ_２］で一様な確率で発生する乱数に応じた値であってもよい。なお、ｗ_１は、０以上ｗ_２以下となり、ｗ_２は、ｗ_１以上１以下となる。

また、ｗは、Ｎ個の変数ｘ_ｉの大きさの平均を表す指標値であってもよい。Ｎ個の変数ｘ_ｉの大きさの平均を表す指標値は、例えば、直前のＮ個の変数ｘ_ｉの二乗平均平方根または絶対値平均である。また、例えば、ｗは、Ｎ個の変数ｘ_ｉの大きさの平均を表す指標値以上且つ１以下の、乱数により定まる値であってもよい。

また、例えば、ｗは、０から１以下の間で初期時刻から終了時刻まで時間経過に従って増加する増加係数であってもよい。増加係数は、例えば、初期時刻で０、終了時刻で１となる時刻を変数とする一次関数、または、一次関数の平方根である。また、例えば、ｗは、増加係数以上、１以下の、乱数により定まる値であってもよい。

さらに、更新によってｘ_ｉ＞１となった場合、変数ｘ_ｉに対応する変数ｙ_ｉに係数ｒを乗算してもよい。すなわち、更新部５０は、値が第１値より小さい第１ベクトルの要素に対応する第２ベクトルの要素、または、第２値より大きい第１ベクトルの要素に対応する第２ベクトルの要素を、もとの前記第２ベクトルの要素に、第２係数を乗じた値に更新するように構成されていてもよい。例えば、更新部５０は、値が－１より小さい第１ベクトルの要素に対応する第２ベクトルの要素、または、値が１より大きい第１ベクトルの要素に対応する第２ベクトルの要素を、もとの第２ベクトルの要素に第２係数を乗じた値に更新するように構成されていてもよい。ここで、第２係数は上述の係数ｒに相当する。

また、更新によって変数ｘ_ｉの絶対値が１より大きくなった場合、変数ｙ_ｉは、０または予め定められた値に変更されてもよい。また、更新によってｘ_ｉ＞１となった場合、変数ｘ_ｉに対応する変数ｙ_ｉの値を擬似乱数に設定してもよい。例えば、［－０．１，０．１］の範囲の乱数を使うことができる。すなわち、更新部５０は、値が第１値より小さい第１ベクトルの要素に対応する第２ベクトルの要素の値、または、値が第２値より大きい第１ベクトルの要素に対応する第２ベクトルの要素の値を、擬似乱数に設定するように構成されていてもよい。

（第１アルゴリズム）
以上のようにして｜ｘ_ｉ｜＞１とならないように更新すると、式（８）の非線形項Ｋｘ^２ _ｉ（ｔ＋Δｔ）を除去しても、ｘ_ｉの値が発散することはなくなる。したがって、式（８）のアルゴリズムに代わり、下記の式（９）の第１アルゴリズムを使うことが可能となる。

なお、ＱＵＢＯ問題の場合、式（９）のｚ_ｉ（ｔ＋Δｔ）は、下記の式（１０）により表される。

上述の式（９）の第１アルゴリズムでは、必ずしも擬似乱数を使う必要がない。また、式（９）の第１アルゴリズムは、式（８）と同様、ハミルトン方程式を解くものであり、変数ｙ_ｉは運動量に相当する。そのため、シンプレクティック・オイラー法を使うことにより、時間ステップΔｔとして小さな値を使わなくても、安定的に解を求めることができる。また、式（９）の第１アルゴリズムにおいても、３次以上の目的関数を有する組合せ最適化問題を解くことが可能である。

図９Ａおよび図９Ｂは、式（９）の第１アルゴリズムを使ってＧ－ｓｅｔのＧ２２を１０００回解いたときの結果を示している。図９Ａおよび図９Ｂでは、式（９）の第１アルゴリズムを使っている。第２係数ｒの値は０に設定され、時間ステップはΔｔ＝１に設定されている。ｗは、区間［ｘ_ａｖ２，１］における一様乱数とした。ｘ_ａｖ２は、Ｎ個の変数ｘ_ｉの二乗平均平方根である。その他の計算条件は、図６Ａおよび図６Ｂと同様であるものとする。なお、式（９）では、非線形項がないため、時間ステップΔｔを図６Ａおよび図６Ｂの倍に設定することができた。

図９Ａは、カット数の平均値を示している。一方、図９Ｂは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Ｃｍａｘの定義は、図６Ａおよび図６Ｂと同様である。また、図９Ａおよび図９Ｂの両グラフにおいて実線で示されたデータは、式（９）の第１アルゴリズムの適用時の結果に相当している。一方、図９Ａおよび図９Ｂの両グラフにおいて破線で示されたデータは、式（８）のアルゴリズムを使ったときの結果に相当する。

図９Ａおよび図９Ｂを参照すると、図６Ａおよび図６Ｂと比べ、カット数の平均値およびカット数の最大値のいずれも最適解に近づいていることがわかる。ただし、図９Ａおよび図９Ｂの結果でも依然として計算値と最適解の間に差がある。この誤差は、第３ベクトルの要素の値ｚ_ｉを定義する基本式において、変数としてスピンｓ_ｉではなく、連続変数ｘ_ｉが使われていることに起因していると考えられる。特に、高次の項が増えるほど、ｚ_ｉにおける変数ｘどうしの積演算は、誤差を大きくする原因となりうる。例えば、１より大きい変数を複数回乗算すると、値が１より著しく大きくなってしまう。

（第２アルゴリズム）
誤差を軽減するために、式（９）の第１アルゴリズムをさらに改良した。具体的には、下記の式（１１）のように、ｚ_ｉにおいて、連続変数ｘ_ｉに代わって連続変数ｘ_ｉを符号関数で変換した値ｓｇｎ（ｘ_ｉ）を代入した。連続変数ｘ_ｉを符号関数で変換した値ｓｇｎ（ｘ_ｉ）は、スピンｓ_ｉに相当する。

なお、ＱＵＢＯ問題の場合、式（１１）のｚ_ｉ（ｔ＋Δｔ）は、下記の式（１２）により表される。

式（１１）では、ｚ_ｉの１階のテンソルを含む項の係数αを定数（例えば、α＝１）にしてもよい。式（１１）の第２アルゴリズムは、式（８）、（９）とは異なり、ハミルトン方程式を解くものではない。式（１１）は、外場によって制御された力学系であると見なすことが可能である。式（１１）の第２アルゴリズムでは、高次の目的関数を有するＨＯＭＯを扱った場合、ｚ_ｉのどのスピンどうしの積も－１または１のいずれかの値をとるため、積演算による誤差の発生を防ぐことができる。

上述の式（１１）の第２アルゴリズムのように、計算サーバ３が計算するデータは、さらに、ｓ_ｉ（ｉ＝１、２、・・・、Ｎ）を要素とする第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）を含んでいてもよい。第１ベクトルのそれぞれの要素を符号関数で変換することにより、第４ベクトルを得ることができる。すなわち、作用演算部５１は、イジングモデルのエネルギー式を、すべての項に含まれる変数について偏微分した形式の基本式を使って、第３ベクトルのそれぞれの要素の値を更新するように構成されていてもよい。ここで、基本式の変数として、第１ベクトルの要素または、第１ベクトルの要素を符号関数で変換した第４ベクトルの要素を使うことができる。

図１０Ａおよび図１０Ｂは、式（１１）の第２アルゴリズムを使ってＧ－ｓｅｔのＧ２２を１０００回解いたときの結果を示している。使用されるアルゴリズムの違いを除けば、図１０Ａおよび図１０Ｂの計算条件（例えば、時間ステップΔｔ、各係数、ｘ_ｉを－１～１の間に保つために定義されるｗ）は、図９Ａおよび図９Ｂと同様であるものとする。図１０Ａは、カット数の平均値を示している。一方、図１０Ｂは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Ｃｍａｘの定義は、図６Ａおよび図６Ｂ、図９Ａおよび図９Ｂと同様である。また、図１０Ａおよび図１０Ｂの両グラフにおいて実線で示されたデータは、式（１１）の第２アルゴリズムの適用時の結果に相当する。一方、図１０Ａおよび図１０Ｂの両グラフにおいて破線で示されたデータは、式（８）のアルゴリズムを使ったときの結果に相当する。

図１０Ａおよび図１０Ｂを参照すると、カット数の平均値およびカット数の最大値のいずれにおいても、図９Ａおよび図９Ｂと比べて最適解により近い値が得られることがわかる。図１０Ｂを参照すると、式（１１）の第２アルゴリズムを使うことにより、カット数の最大値１３３５９が得られることがわかる。

（第３アルゴリズム）
式（９）の第１アルゴリズムをさらに下記の式（１３）のように変形してもよい。

なお、ＱＵＢＯ問題の場合、式（１３）のｚ_ｉ（ｔ＋Δｔ）は、下記の式（１４）により表される。

式（１３）の第３アルゴリズムは、外力に相当する項ｆ_ｉの計算方法が上述の各例とは異なっている。（１３）４番目の式を使って計算された値ｚ_ｉを符号関数で変換し、１で規格化している。すなわち、作用演算部５１は、第１ベクトルの要素を変数として計算した基本式の値（ｚ_ｉ）を第１関数で変換した値に基づいて第３ベクトルのそれぞれの要素を更新するように構成されていてもよい。第１関数として、例えば、符号関数を使うことができる。ただし、後述するように、その他の関数を第１関数として使ってもよい。

また、式（１３）では、係数ｃの代わりに、関数ｇ（ｔ）が使われている。一般に、第３ベクトルの要素の値ｚ_ｉの計算結果への寄与度は、問題によって異なる。しかし、式（１３）では、第３ベクトルの要素の値ｚ_ｉが１で規格化されているため、問題ごとに係数ｃの値を決定する必要がなくなる。関数ｇ（ｔ）として、例えば、下記の式（１５）を使うことができる。

式（１５）の関数は、更新回数に応じて、単調増加してから単調減少する。ただし、上述の式（１５）は一例にすぎず、ｇ（ｔ）として、ｐ（ｔ）をパラメータとする、これとは異なる関数を使ってもよい。すなわち、作用演算部５１は、第１係数ｐをパラメータとする第２関数を乗じることによって第３ベクトルのそれぞれの要素を更新するように構成されていてもよい。

図１１Ａおよび図１１Ｂは、式（１３）の第３アルゴリズムを使ってＧ－ｓｅｔのＧ２２を１０００回解いたときの結果を示している。使用されるアルゴリズムの違いを除けば、図１１Ａおよび図１１Ｂの計算条件（例えば、時間ステップΔｔ、使われる係数、ｘ_ｉを－１～１の間に保つために定義されるｗ）は、図９Ａおよび図９Ｂと同様であるものとする。図１１Ａは、カット数の平均値を示している。一方、図１１Ｂは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Ｃｍａｘの定義は、図６Ａおよび図６Ｂ、図９Ａおよび図９Ｂと同様である。また、図１１Ａおよび図１１Ｂの両グラフにおいて実線で示されたデータは、式（１３）の第３アルゴリズムの適用時の結果に相当する。一方、図１１Ａおよび図１１Ｂの両グラフにおいて破線で示されたデータは、式（８）のアルゴリズムを使ったときの結果に相当する。

図１１Ａおよび図１１Ｂを参照すると、カット数の平均値およびカット数の最大値のいずれにおいても、式（８）のアルゴリズムと比べて最適解により近い値が得られることがわかる。図１１Ｂを参照すると、式（１３）の第３アルゴリズムを使うことにより、最大値１３３５９が得られていることがわかる。

（変形例）
式（９）、式（１１）および式（１３）のアルゴリズムでは、基本式（ｚ_ｉの式）の１階のテンソルを含む項の係数αを定数係数（例えば、α＝１）として計算を行ってもよい。また、式（９）、式（１１）および式（１３）のアルゴリズムでは、基本式（ｚ_ｉの式）の１階のテンソルを含む項の係数αとして更新回数に応じて単調減少または単調増加する係数を使ってもよい。この場合、基本式の１階のテンソルを含む項は、更新回数に応じて単調減少または単調増加する。

上述の式（９）の第１アルゴリズムおよび式（１１）の第２アルゴリズムには、係数ｃが含まれている。係数ｃをＪ^（２）行列の最大固有値の逆数に近い値にしたい場合、Ｊ^（２）行列の最大固有値を計算するか、Ｊ^（２）行列の最大固有値の見積もりを行う必要がある。最大固有値の計算は、必要な計算量が大きい。一方、最大固有値の見積もりは、値の正確性が保証されていない。そこで、係数ｃの代わりに上述の式（１５）のような、更新回数に応じて値が変動する関数を使うことができる。また、係数ｃの代わりに、下記の式（１６）のように、第１のベクトル（ｘ_１、ｘ_２、・・・、ｘ_Ｎ）および第３ベクトル（ｚ_１、ｚ、・・・、ｚ_Ｎ）に基づいて計算される、近似値ｃ１を使ってもよい。

式（１６）を参照すると、分母と分子がいずれもベクトルのノルムとなっている。式（１６）のように、ベクトルのノルムとして、ベクトルの各要素の２乗和の平方根である、Ｌ２ノルムを使うことができる。ただし、ベクトルの要素の絶対値の和である、Ｌ１ノルムなど、その他の定義によるノルムを使ってもよい。

すなわち、更新部５０は、第１ベクトルのノルムを第３ベクトルのノルムで除算した第３係数ｃ１を計算し、第２ベクトルの要素に、第１係数ｐ（ｔ＋Δｔ）と更新された第１ベクトルの対応する要素との積を重み付けした値ならびに、第３ベクトルの対応する要素を第３係数ｃ１で重み付けした値を加算することによって第２ベクトルの要素を更新するように構成されていてもよい。

さらに、係数ｃの代わりに、下記の式（１７）のような内積によって定義される近似値ｃ´１を使ってもよい。

すなわち、更新部５０は、第１ベクトルどうしの内積を、第１ベクトルと第３ベクトルの内積の絶対値で除算した第３係数ｃ´１を計算し、第２ベクトルの要素に、第１係数ｐ（ｔ＋Δｔ）と更新された第１ベクトルの対応する要素との積を重み付けした値ならびに、第３ベクトルの対応する要素を第３係数ｃ´１で重み付けした値を加算することによって第２ベクトルの要素を更新するように構成されていてもよい。

近似値ｃ１、ｃ´１は、各計算タイミングにおける第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）および第３ベクトル（ｚ_１，ｚ_２，・・・，ｚ_Ｎ）の値に基づいて計算されるため、係数ｃのような定数ではなく、動的に制御される係数となる。なお、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）および第３ベクトル（ｚ_１，ｚ_２，・・・，ｚ_Ｎ）については、変数の更新処理で計算されたものを利用することができるため、近似値ｃ１、ｃ´１を計算したとしても、計算量が大幅に増えることはない。局所磁場のないイジング問題において、（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）がＪ^（２）の最大固有値に対応する固有ベクトルであるとき、近似値ｃ１、ｃ´１はＪ^（２）の最大固有値の逆数に等しくなる。また、（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）が固有ベクトルからずれていると、近似値ｃ１、ｃ´１はＪ^（２）の最大固有値の逆数より大きな値となるため、解の収束が早まる。

図１２Ａおよび図１２Ｂは、式（９）の第１アルゴリズムにおいて、係数ｃの代わりに近似値ｃ１を使い、Ｇ－ｓｅｔのＧ２２を１０００回解いたときの結果を示している。図１３Ａおよび図１３Ｂは、式（１１）の第２アルゴリズムにおいて、係数ｃの代わりに近似値ｃ１を使い、Ｇ－ｓｅｔのＧ２２を１０００回解いたときの結果を示している。

図１２Ａおよび図１３Ａは、カット数の平均値を示している。一方、図１２Ｂおよび図１３Ｂは、カット数の最大値を示している。各軸の対応関係および両グラフにある水平方向の破線Ｃｍａｘの定義は、上述の各グラフと同様である。図１２Ａ、図１２Ｂ、図１３Ａおよび図１３Ｂにおいて破線で示されたデータは、式（８）のアルゴリズムを使ったときの結果を示している。

図１２Ａ、図１２Ｂ、図１３Ａおよび図１３Ｂを参照すると、カット数の平均値およびカット数の最大値のいずれにおいても、式（８）のアルゴリズムと比べて最適解により近い値が得られることがわかる。特に、式（１１）の第２アルゴリズムにおいて、カット数の最大値１３３５９が得られていることがわかる。

なお、式（９）の第１アルゴリズムおよび式（１１）の第２アルゴリズムでは、近似値ｃ１およびｃ´１の代わりに下記の式（１８）で定義される近似値ｃ２またはｃ´２を使ってもよい。

すなわち、更新部５０は、第１ベクトルのそれぞれの要素を符号関数で変換した第４ベクトルのノルムを、第３ベクトルのノルムで除算した第３係数ｃ２を計算し、第２ベクトルの要素に、第１係数ｐ（ｔ＋Δｔ）と更新された第１ベクトルの対応する要素との積を重み付けした値ならびに、第３ベクトルの対応する要素を第３係数ｃ２で重み付けした値を加算することによって第２ベクトルの要素を更新するように構成されていてもよい。

また、更新部５０は、第１ベクトルのそれぞれの要素を符号関数で変換した第４ベクトルどうしの内積を、第４ベクトルと第３ベクトルの内積の絶対値で除算した第３係数ｃ´２を計算し、第２ベクトルの要素に、第１係数ｐ（ｔ＋Δｔ）と更新された第１ベクトルの対応する要素との積を重み付けした値ならびに、第３ベクトルの対応する要素を第３係数ｃ´２で重み付けした値を加算することによって第２ベクトルの要素を更新するように構成されていてもよい。

式（１８）の第３ベクトル（ｚ_１，ｚ_２，・・・，ｚ_Ｎ）については、アルゴリズムで計算されたものを使うことができるため、近似値ｃ２、ｃ´２を求めたとしても、計算量が大幅に増えることはない。

アルゴリズムの実行中における各ベクトルの値を使って、近似値ｃ１、ｃ´１、ｃ２、ｃ´２を計算すると、各計算タイミングによって値が激しく変動することがある。近似値ｃ１、ｃ´１、ｃ２、ｃ´２の変動を抑制するため、近似値ｃ１、ｃ´１、ｃ２、ｃ´２の代わりに近似値ｃ１、ｃ´１、ｃ２、ｃ´２を所定の規則に基づいて変換した値を使ってもよい。例えば、所定の規則として下記の式（１９）を使うことができる。

ここで、γに１より小さい値を設定することができる。式（１９）のｃ（ｔ＋Δｔ）には、例えば、上述の式（１６）～（１８）によって計算された近似値が代入される。ｃ（ｔ＋Δｔ）を、各計算タイミングにおいて振動成分を含む信号をサンプリングした値であるとみなすと、ｄ（ｔ＋Δｔ）は、ｃ（ｔ＋Δｔ）が一定帯域のローパスフィルタを通過した後の値に相当するといえる。

すなわち、更新部５０は、第３係数（近似値ｃ１、ｃ´１、ｃ２、ｃ´２のいずれか）がローパスフィルタを通過した後の値である第４係数を計算し、第３係数に代わり、第４係数を用いて第２ベクトルの要素を更新するように構成されていてもよい。

上述では、シミュレーテッド分岐アルゴリズムを用いて、イジングモデルの解を求める例について説明した。ただし、シミュレーテッド分岐アルゴリズムによる求解が可能な組合せ最適化問題は、イジング問題に限られない。シミュレーテッド分岐アルゴリズムを用いることによって、一般的な２値変数の組合せ最適化問題を解くことが可能である。例えば、上述の各アルゴリズムは、目的関数の変数が、ａ（第１値）と、ａより大きいｂ（第２値）のいずれかをとる２値変数である、組合せ最適化問題に適用することが可能である。また、一定の更新回数の後に目的関数の解を求める場合、符号関数の代わりに、値域がａまたはｂの２値である関数ｆ（ｘ_ｉ）を使ってもよい。この関数ｆ（ｘ_ｉ）がとる値は、変数ｘ_ｉの値をしきい値ｖ（ａ＜ｖ＜ｂ）と比較した結果に基づいて決まる。例えば、ｘ_ｉ＜ｖであるならば、ｆ（ｘ_ｉ）＝ａとなる。また、ｖ＜ｘ_ｉであるならば、ｆ（ｘ_ｉ）＝ｂとなる。例えば、ｘ_ｉ＝ｖである場合、ｆ（ｘ_ｉ）＝ａまたは、ｆ（ｘ_ｉ）＝ｂとなる。ここで、しきい値ｖの値として、例えば、（ａ＋ｂ）／２を使うことができる。上述の関数ｆ（ｘ_ｉ）は、第１ベクトルの要素を第４ベクトルの要素に変換する第１関数として使われてもよい。

また、例えば、上述の式（９）の第１アルゴリズム、式（１１）の第２アルゴリズムおよび式（１３）の第３アルゴリズムを使った場合、更新によって変数ｘ_ｉがａより小さくなったとき、変数ｘ_ｉの値を｛ｖ－ｗ_－｝に変更する。ｗ_－は、０以上（ｖ－ａ）以下の実数である。また、更新によって変数ｘ_ｉがｂより大きくなったとき、変数ｘ_ｉの値を｛ｖ＋ｗ_＋｝に変更する。ｗ_＋は、０以上（ｂ－ｖ）以下の実数である。

例えば、ｗ_－は、０以上（ｖ－ａ）以下の予め定められた値であってもよい。ｗ_＋は、０以上（ｂ－ｖ）以下の予め定められた値であってもよい。

また、ｗ_－は、０以上（ｖ－ａ）以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値であってもよい。また、ｗ_＋は、０以上（ｂ－ｖ）以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値であってもよい。

また、ｗ_－およびｗ_＋は、複数の第１変数（ｘ_ｉ）のしきい値（ｖ）からの偏差の大きさに対する平均を表す指標値（ｘ_ａｖｅ）であってもよい。例えば、指標値（ｘ_ａｖｅ）は、複数の要素の第１変数（ｘ_ｉ）のしきい値（ｖ）からの偏差の二乗平均平方根または絶対値平均である。

また、ｗ_－は、指標値（ｘ_ａｖｅ）以上、（ｖ－ａ）以下の乱数により定まる値であってもよい。また、ｗ_＋は、指標値（ｘ_ａｖｅ）以上、（ｂ－ｖ）以下の乱数により定まる値であってもよい。なお、指標値（ｘ_ａｖｅ）が（ｖ－ａ）を超える場合には、ｗ_－は、（ｖ－ａ）とする。また、指標値が（ｂ－ｖ）を超える場合には、ｗ_＋は、（ｂ－ｖ）とする。

また、ｗ_－およびｗ_＋は、増加係数であってもよい。増加係数は、更新処理の初期時刻において０であり、初期時刻から終了時刻まで時間経過に従って増加する。また、ｗ_－は、増加係数以上、（ｖ－ａ）以下の乱数により定まる値であってもよい。

また、ｗ_＋は、増加係数以上、（ｂ－ｖ）以下の乱数により定まる値であってもよい。なお、増加係数が（ｖ－ａ）を超える場合には、ｗ_－は、（ｖ－ａ）とする。また、増加係数が（ｂ－ｖ）を超える場合には、ｗ_＋は、（ｂ－ｖ）とする。

以上から、目的関数の変数が第１値（ａ）と第２値（ｂ）のいずれかをとる離散変数である組合せ最適化問題に対して、上述の式（９）の第１アルゴリズム、式（１１）の第２アルゴリズムおよび式（１３）の第３アルゴリズムを適用した場合、第１変数ｘ_ｉおよび第２変数ｙ_ｉの更新処理において、更新部５０は、次のような処理を行う。

すなわち、更新部５０は、第１変数（ｘ_ｉ）が第１値（ａ）より小さい場合、第１変数（ｘ_ｉ）を、第１値（ａ）以上、しきい値（ｖ）以下の値に変更する。さらに、更新部５０は、第１変数（ｘ_ｉ）が第２値（ｂ）より大きい場合、第１変数（ｘ_ｉ）を、予め定められたしきい値（ｖ）以上、第２値（ｂ）以下の値に変更する。

より具体的には、例えば、更新部５０は、第１変数（ｘ_ｉ）が第１値（ａ）より小さい場合、第１変数（ｘ_ｉ）を、第１値（ａ）以上しきい値（ｖ）以下の予め定められた値、または、第１値（ａ）以上しきい値（ｖ）以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ）が第２値（ｂ）より大きい場合、第１変数（ｘ_ｉ）を、しきい値（ｖ）以上第２値（ｂ）以下の予め定められた値、または、しきい値（ｖ）以上第２値（ｂ）以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値に変更してもよい。

また、例えば、更新部５０は、第１変数（ｘ_ｉ）が第１値（ａ）より小さい場合、第１変数（ｘ_ｉ）を、しきい値（ｖ）から指標値を減算した値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ）が第２値（ｂ）より大きい場合、第１変数（ｘ_ｉ）を、しきい値（ｖ）に指標値を加算した値に変更してもよい。この場合、指標値は、複数の要素の第１変数（ｘ_ｉ）のしきい値（ｖ）からの偏差の大きさに対する平均を表す。例えば、指標値は、複数の要素の第１変数（ｘ_ｉ）のしきい値（ｖ）からの偏差の二乗平均平方根または絶対値平均である。

また、例えば、更新部５０は、第１変数（ｘ_ｉ）が第１値（ａ）より小さい場合、第１変数（ｘ_ｉ）を、第１値（ａ）以上、しきい値（ｖ）から指標値を減じた値以上の、乱数により定まる値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ）が第２値（ｂ）より大きい場合、第１変数（ｘ_ｉ）を、しきい値（ｖ）に指標値を加算した値以上、第２値（ｂ）以下の、乱数により定まる値に変更してもよい。

また、更新部５０は、第１変数（ｘ_ｉ）が第１値（ａ）より小さい場合、第１変数（ｘ_ｉ）を、しきい値（ｖ）から増加係数を減算した値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ）が第２値（ｂ）より大きい場合、第１変数（ｘ_ｉ）を、しきい値（ｖ）に増加係数を加算した値に変更してもよい。この場合、増加係数は、初期時刻において０であり、初期時刻から終了時刻まで時間経過に従って増加する。

また、更新部５０は、第１変数（ｘ_ｉ）が第１値（ａ）より小さい場合、第１変数（ｘ_ｉ）を、第１値（ａ）以上、しきい値（ｖ）に増加係数を加算した値以下における、乱数により定まる値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ）が第２値（ｂ）より大きい場合、第１変数（ｘ_ｉ）を、しきい値（ｖ）に増加係数を加算した値以上、第２値（ｂ）以下における、乱数により定まる値に変更してもよい。

ここまでは、シンプレクティック・オイラー法によって実装されたシミュレーテッド分岐アルゴリズムの例と、それぞれのアルゴリズムを使って組合せ最適化問題を計算した結果について説明した。以降では、上述のアルゴリズムの実装例について述べる。

（ＰＣクラスタへの実装例）
はじめに、ＰＣクラスタへ上述のアルゴリズムを実装した例について説明する。ＰＣクラスタとは、複数台のコンピュータを接続し、１台のコンピュータでは得られない計算性能を実現するシステムである。例えば、図１に示した情報処理システム１００は、複数台の計算サーバ３およびプロセッサを含んでおり、ＰＣクラスタとして利用することが可能である。例えば、ＰＣクラスタにおいては、ＭＰＩ（Ｍｅｓｓａｇｅ Ｐａｓｓｉｎｇ Ｉｎｔｅｒｆａｃｅ）を使うことにより、情報処理システム１００のような複数の計算サーバ３にメモリが分散して配置されている構成でも並列的な計算を実行することが可能である。例えば、ＭＰＩを使って管理サーバ１の制御プログラム１４Ｅ、各計算サーバ３の計算プログラム３４Ｂおよび制御プログラム３４Ｃを実装することができる。

ＰＣクラスタで利用するプロセッサ数がＱである場合、それぞれのプロセッサに、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）に含まれる変数ｘ_ｉのうち、Ｌ個の変数の計算を行わせることができる。同様に、それぞれのプロセッサに、第２ベクトル（ｙ_１，ｙ_２，・・・，ｙ_Ｎ）に含まれる変数ｙ_ｉのうち、Ｌ個の変数の計算を行わせることができる。すなわち、プロセッサ＃ｊ（ｊ＝１，２，・・・，Ｑ）は、変数｛ｘ_ｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝および｛ｙ_ｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の計算を行う。また、プロセッサ＃ｊによる｛ｙ_ｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の計算に必要な下記の式（２０）に示されたテンソルＪ^（ｎ）は、プロセッサ＃ｊがアクセス可能な記憶領域（例えば、レジスタ、キャッシュ、メモリなど）に保存されるものとする。

ここでは、それぞれのプロセッサが第１ベクトルおよび第２ベクトルの一定数の変数を計算する場合を説明した。ただし、プロセッサによって、計算する第１ベクトルおよび第２ベクトルの変数の数が異なっていてもよい。例えば、計算サーバ３に実装されるプロセッサによって性能差がある場合、プロセッサの性能に応じて計算対象とする変数の数を決めることができる。

すなわち、情報処理装置（例えば、計算サーバ３）は、複数のプロセッサを備えていてもよい。更新部５０は、複数のプロセッサを含んでおり、更新部５０の複数のプロセッサのそれぞれは、第１ベクトルの一部の要素の値および第２ベクトルの一部の要素の値を更新するように構成されていてもよい。

変数ｙ_ｉの値を更新するためには、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）または第１ベクトルの各要素を２値変数に変換した第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）のすべての成分の値が必要となる。２値変数への変換は、例えば、符号関数ｓｇｎ（）を使うことによって行うことができる。そこで、Ａｌｌｇａｔｈｅｒ関数を使い、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）または第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）のすべての成分の値をＱ個のプロセッサに共有させることができる。第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）または第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）については、プロセッサ間での値の共有が必要であるものの、第２ベクトル（ｙ_１，ｙ_２，・・・，ｙ_Ｎ）およびテンソルＪ^（ｎ）については、プロセッサ間での値の共有を行うことは必須ではない。プロセッサ間でのデータの共有は、例えば、プロセッサ間通信を使ったり、共有メモリにデータを保存したりすることによって実現することができる。

プロセッサ＃ｊは、変数｛ｚ_ｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の値を計算する。そして、プロセッサ＃ｊは、計算した｛ｚ_ｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の値に基づき、変数｛ｙ_ｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝を更新する。

上述の各式に示したように、ベクトル（ｚ_１，ｚ_２，・・・，ｚ_Ｎ）の計算では、テンソルＪ^（ｎ）と、ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）または（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）との積の計算を含む、積和演算が必要である。積和演算は、上述のアルゴリズムにおいて最も計算量の大きい処理であり、計算速度の向上においてボトルネックとなりうる。そこで、ＰＣクラスタの実装では、積和演算を、Ｑ＝Ｎ／Ｌ個のプロセッサに分散して並列的に実行し、計算時間の短縮をはかることができる。

すなわち、情報処理装置（例えば、計算サーバ３）は、複数のプロセッサを備えていてもよい。作用演算部５１は、複数のプロセッサを含み、作用演算部５１の複数のプロセッサのそれぞれは、第３ベクトルの一部の要素を更新するように構成されていてもよい。更新部５０は、複数のプロセッサを含み、更新部５０の複数のプロセッサのそれぞれは、第１ベクトルの一部の要素および第２ベクトルの一部の要素を更新するように構成されていてもよい。

図１４は、マルチプロセッサ構成の例を概略的に示している。図１４の複数の計算ノードは、例えば、情報処理システム１００の複数の計算サーバ３に相当する。また、図１４の高速リンクは、例えば、情報処理システム１００のケーブル４ａ～４ｃおよびスイッチ５によって形成された計算サーバ３間のインターコネクトに相当する。図１４の共有メモリは、例えば、共有メモリ３２に相当する。図１４のプロセッサは、例えば、各計算サーバ３のプロセッサ３３Ａ～３３Ｄに相当している。なお、図１４には複数の計算ノードが示されているが、単一計算ノードの構成を用いることを妨げるものではない。

図１４には、各構成要素に配置されるデータおよび構成要素間で転送されるデータが示されている。各プロセッサでは、変数ｘ_ｉ、（ｓ_ｉ）、ｙ_ｉ、ｚ_ｉの値が計算される。また、プロセッサと共有メモリ間では、変数ｘ_ｉまたはｓ_ｉが転送される。各計算ノードの共有メモリには、例えば、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）、第２ベクトル（ｙ_１，ｙ_２，・・・，ｙ_Ｎ）のＬ個の変数、およびテンソルＪ^（ｎ）の一部が保存される。なお、式（１１）のアルゴリズムを実行する場合、各計算ノードの共有メモリには、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）に代わって第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）が保存されていてもよい。そして、計算ノード間を接続する高速リンクでは、例えば、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）が転送される。Ａｌｌｇａｔｈｅｒ関数を使った場合、各プロセッサで変数ｙ_ｉおよびｚ_ｉを更新するために、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）の全要素が必要だからである。なお、式（１１）のアルゴリズムにしたがって変数ｚ_ｉを更新する場合には、各プロセッサは第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）の全要素にアクセスする必要がある。このため、高速リンクでは、第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）が転送されてもよい。

ただし、図１４に示したデータの配置および転送は一例にしかすぎない。例えば、各プロセッサが積和演算を含む｛ｚ_ｍ｜ｍ＝（ｊ－１）Ｌ＋１，（ｊ－１）Ｌ＋２，・・・，ｊＬ｝の計算を並列的に実行しているのであれば、それぞれのプロセッサと共有メモリ間および計算ノード間で変数ｚ_ｉの値を転送し、共有されたベクトル（ｚ_１，ｚ_２，・・・，ｚ_Ｎ）を参照し、変数ｙ_ｉの値を計算してもよい。このように、ＰＣクラスタにおけるデータの配置方法、転送方法および並列化の実現方法については、特に問わない。

すなわち、情報処理装置（例えば、計算サーバ３）は、複数のプロセッサからアクセス可能に構成されている共有メモリを備えていてもよい。この場合、更新部５０は、更新された後の第１ベクトルの要素または、更新された後の第１ベクトルのそれぞれの要素を２値変数に変換した第４ベクトルを共有メモリに保存することができる。

次に、ＰＣクラスタに上述の各アルゴリズムを実行させたときにおける結果について説明する。図１５Ａおよび図１５Ｂは、ＰＣクラスタを使ってＮ=３６００の（局所磁場のない）全結合イジング問題を解いたときの結果を示している。図１５Ａおよび図１５Ｂの全結合イジング問題では、結合係数の行列Ｊの各成分の値は、［－１，１］の範囲の一様乱数に設定した。また、合計時間ステップ数は、１００００とした。図１５Ａは、各アルゴリズムを使って１０回全結合イジング問題を解いたときにおけるエネルギーＥ_{Ｉｓｉｎｇ}の平均値を示している。図１５Ｂは、各アルゴリズムを使って１０回全結合イジング問題を解いたときにおける計算時間の平均値を秒単位で示している。

図１５Ａおよび図１５Ｂの棒グラフには、左側から右側に向かって、
（ｉ） 式（８）のアルゴリズムが使われた場合、
（ｉｉ） 式（９）の第１アルゴリズムが使われた場合、
（ｉｉｉ） 式（１１）の第２アルゴリズムが使われた場合、
（ｉｖ） 式（１３）の第３アルゴリズムが使われた場合、
（ｖ） 式（９）の第１アルゴリズムにおいて、係数ｃの代わりに近似値ｃ１が使われた場合、
（ｖｉ） 式（１１）の第２アルゴリズムにおいて、係数ｃの代わりに近似値ｃ１が使われた場合、
の６のケースの結果が示されている。

また、図１５Ａおよび図１５Ｂの棒グラフには、各ケースのそれぞれについて、左側から右側に向かって、プロセッサ数Ｑが１および３６の場合における結果が示されている。式（８）のアルゴリズムが使われた場合、式（９）の第１アルゴリズムが使われた場合、式（１１）の第２アルゴリズムが使われた場合において、係数ｃとして定数０．５Ｄ√（３／Ｎ）を用いた。

図１５Ａの棒グラフを参照すると、（ｉｉ）～（ｖｉ）のケースでは、（ｉ）のケースと比べ、エネルギーＥ_{Ｉｓｉｎｇ}の平均値が低くなっており、最適解により近い解が得られやすいことがわかる。また、図１５Ｂの棒グラフを参照すると、マルチプロセッサ構成による並列計算を行うことにより、計算時間が著しく短縮されることがわかる。

（ＧＰＵへの実装例）
ＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）を使って上述の各アルゴリズムの計算を行ってもよい。図１６は、ＧＰＵを使った構成の例を概略的に示している。図１６には、互いに高速リンクで接続された複数のＧＰＵが示されている。それぞれのＧＰＵには、共有メモリにアクセス可能な複数のコアが搭載されている。また、図１６の構成例では、複数のＧＰＵが高速リンクを介して接続されており、ＧＰＵクラスタを形成している。例えば、ＧＰＵが図１のそれぞれの計算サーバ３に搭載されている場合、高速リンクは、ケーブル４（４ａ～４ｃ）およびスイッチ５によって形成された計算サーバ３間のインターコネクトに相当する。なお、図１６の構成例では、複数のＧＰＵが使われているが、ひとつのＧＰＵを使った場合にも、並列的な計算を実行することが可能である。すなわち、図１６のそれぞれのＧＰＵは、図１４のそれぞれの計算ノードに相当する計算を実行できる。すなわち、情報処理装置（計算サーバ３）のプロセッサは、Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ（ＧＰＵ）のコアであってもよい。

ＧＰＵにおいて、変数ｘ_ｉおよびｙ_ｉ、ならびにテンソルＪ^（ｎ）はデバイス変数として定義される。ＧＰＵは、変数ｙ_ｉの更新に必要なテンソルＪ^（ｎ）と第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）または第４ベクトル（ｓ_１，ｓ_２，・・・，ｓ_Ｎ）の積を、行列ベクトル積関数によって並列的に計算することができる。なお、行列とベクトルの積演算を繰り返し実行することにより、テンソルとベクトルの積を求めることができる。また、第１ベクトル（ｘ_１，ｘ_２，・・・，ｘ_Ｎ）の計算と、第２ベクトル（ｙ_１，ｙ_２，・・・，ｙ_Ｎ）のうち、積和演算以外の部分については、それぞれのスレッドにｉ番目の要素（ｘ_ｉ，ｙ_ｉ）の更新処理を実行させ、処理の並列化を実現することができる。

図１７Ａおよび図１７Ｂは、ＧＰＵを使ってＮ=３６００の全結合イジング問題を解いたときの結果を示している。図１７Ａの棒グラフは、各アルゴリズムを使って１０回全結合イジング問題を解いたときにおけるエネルギーＥ_{Ｉｓｉｎｇ}の平均値を示している。図１７Ｂの棒グラフは、各アルゴリズムを使って１０回全結合イジング問題を解いたときにおける計算時間の平均値を秒単位で示している。

図１７Ａおよび図１７Ｂの棒グラフには、左側から右側に向かって、式（８）のアルゴリズムが使われた場合、式（９）の第１アルゴリズムが使われた場合、式（１１）の第２アルゴリズムが使われた場合、式（１３）の第３アルゴリズムが使われた場合の４つのケースにおける結果が示されている。式（８）のアルゴリズムが使われた場合、式（９）の第１アルゴリズムが使われた場合、式（１１）の第２アルゴリズムが使われた場合において、係数ｃとして定数０．５Ｄ√（３／Ｎ）を用いた。いずれの結果においても合計時間ステップ数は、１００００である。また、それぞれのアルゴリズムについて、左側が１計算ノードのＰＣクラスタ、右側が１ＧＰＵを用いたときの結果が示されている。

図１７Ａの棒グラフを参照すると、式（９）、式（１１）および式（１３）のアルゴリズムが使われた場合には、式（８）のアルゴリズムが使われた場合と比べ、エネルギーＥ_{Ｉｓｉｎｇ}の平均値が低くなっており、最適解により近い解が得られやすいことがわかる。また、図１７Ｂの棒グラフを参照すると、ＧＰＵによる並列計算を行うことにより、１計算ノードのＰＣクラスタと比べて計算時間を著しく短縮できることがわかる。これは、ＧＰＵの計算の並列度が、一般的なＣＰＵと比べて高いことに起因する。

（第１アルゴリズムの処理フローの第１例）
図１８は、第１アルゴリズムを実行する場合における情報処理システム１００の処理の流れの第１例を示す。情報処理システム１００は、式（９）および式（１０）に示す第１アルゴリズムを用いて最適化問題を解く場合、例えば、図１８に示す流れで処理を実行する。

まず、Ｓ１０１において、更新部５０は、パラメータを設定する。具体的には、更新部５０は、Ｎ×Ｎ個の結合係数を含む行列であるＪ、および、Ｎ個の局所磁場を表す局所磁場係数を含む配列であるｈを設定する。なお、更新部５０は、ＨＯＢＯ問題を解く場合には、Ｊおよびｈに代えて、Ｎ^ｎ個の作用係数を含むｎ次テンソルであるＪ^（ｎ）を設定する。この場合、ｎは、ＨＯＢＯ問題の目的関数の変数の次数を表す。さらに、更新部５０は、係数であるＤ、係数であるｃ、単位時間を表すΔｔ、終了時刻を表すＴ、関数であるｐ（ｔ）、および、関数であるα（ｔ）を設定する。ｐ（ｔ）およびα（ｔ）は、ｔ＝初期時刻（例えば０）で０、ｔ＝終了時刻（Ｔ）で１となる増加関数である。更新部５０は、Ｊ、ｈをユーザからの受け取った情報に応じて設定する。更新部５０は、Ｄ、ｃ、Δｔ、Ｔ、ｐ（ｔ）およびα（ｔ）を、ユーザから受け取ったパラメータに応じて設定してもよいし、予め決定されており変更できないパラメータを設定してもよい。

続いて、Ｓ１０２において、更新部５０は、変数を初期化する。具体的には、更新部５０は、時刻を表す変数であるｔを初期時刻（例えば、０）に初期化する。さらに、更新部５０は、Ｎ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））のそれぞれおよびＮ個の第２変数（ｙ_１～ｙ_Ｎ）のそれぞれに、ユーザから受け取った初期値、予め定められた固定値、または、乱数を代入する。

続いて、更新部５０は、Ｓ１０３とＳ１１８との間のループ処理を、ｔがＴより大きくなるまで繰り返す。１回のループ処理において、更新部５０は、対象時刻（ｔ＋Δｔ）におけるＮ個の第１変数（ｘ_１（ｔ＋Δｔ）～ｘ_Ｎ（ｔ＋Δｔ））を、直前時刻（ｔ）におけるＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））に基づき算出する。また、１回のループ処理において、更新部５０は、対象時刻（ｔ＋Δｔ）におけるＮ個の第２変数（ｙ_１（ｔ＋Δｔ）～ｙ_Ｎ（ｔ＋Δｔ））を、直前時刻（ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））に基づき算出する。

なお、直前時刻（ｔ）は、対象時刻（ｔ＋Δｔ）より単位時間（Δｔ）前の時刻である。すなわち、更新部５０は、Ｓ１０３とＳ１１８との間のループ処理を繰り返すことにより、Ｎ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎に順次に更新する。

Ｓ１０４において、更新部５０は、直前時刻（ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））の大きさの平均を表す指標値（ｘ_ａｖｅ）を算出する。例えば、指標値（ｘ_ａｖｅ）は、直前時刻におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））の二乗平均平方根または絶対値平均である。例えば、更新部５０は、二乗平均平方根を算出する場合、式（２１－１）に示す演算を実行する。また、例えば、更新部５０は、絶対値平均を算出する場合、式（２１－２）に示す演算を実行する。なお、更新部５０は、後述するＳ１０９において、指標値（ｘ_ａｖｅ）を用いない場合、Ｓ１０４の処理を実行しない。

続いて、更新部５０は、Ｓ１０５とＳ１１１との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。ｉは、１からＮまでの整数であり、Ｎ個の要素のうちの処理対象を表すインデックスである。Ｎ個の要素のそれぞれは、第１変数（ｘ_ｉ（ｔ））および第２変数（ｙ_ｉ（ｔ））が対応付けられる。Ｓ１０５とＳ１１１との間のループ処理において、更新部５０は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。

Ｓ１０６において、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））に、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））と予め定められた定数（Ｄ）と単位時間（Δｔ）とを乗算した値を加算することにより算出する。具体的には、更新部５０は、式（２２）を算出する。

続いて、Ｓ１０７において、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が、予め定められた第２値（＋１）より大きいか否かを判断する。本例において、第２値は、＋１である。第２値は、連続量である第１変数（ｘ_ｉ（ｔ））の単位量である。更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が第２値以下である場合、処理をＳ１１１に進める（Ｓ１０７のＮｏ）。更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が第２値より大きい場合、処理をＳ１０８に進める。

続いて、Ｓ１０８において、更新部５０は、０以上、第２値（例えば＋１）以下の範囲内の所定区間で一様な確率で発生する乱数（ｒ）を生成する。なお、更新部５０は、後述するＳ１０９において、乱数（ｒ）を用いない場合、Ｓ１０８の処理を実行しない。

続いて、Ｓ１０９において、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））の制約処理をする。具体的には、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、符号を変更せずに、絶対値を０以上、第２値以下の値に変更する。本例においては、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））が第１値（ａ）である－１より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、－１以上、しきい値（ｖ）である０以下の値に変更する。さらに、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が第２値（ｂ）である＋１より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０以上、＋１以下の値に変更する。

例えば、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が－１より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、－１以上０以下の予め定められた値、または、Ｓ１０８で生成した乱数（ｒ）に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が＋１より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０以上＋１以下の予め定められた値、または、０からＳ１０８で生成した乱数（ｒ）を減じた値に変更してもよい。

また、例えば、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が－１より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０からＳ１０４で算出した指標値（ｘ_ａｖｅ）を減算した値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が＋１より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０に指標値（ｘ_ａｖｅ）を加算した値に変更してもよい。

また、例えば、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が－１より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、－１以上、０から指標値（ｘ_ａｖｅ）を減じた値以上の、Ｓ１０８で算出した乱数（ｒ）により定まる値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が＋１より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０に指標値（ｘ_ａｖｅ）を加算した値以上、＋１以下の、Ｓ１０８で算出した乱数（ｒ）により定まる値に変更してもよい。この場合、更新部は、例えば、式（２３）を算出する。

また、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が－１より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０から増加係数を減算した値に変更してもよい。増加関数は、０から＋１以下の初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで時間経過に従って増加する増加係数に変更してもよい。この場合、更新部５０は、初期時刻（ｔ＝０）で０、終了時刻（ｔ＝Ｔ）で第２値（例えば＋１）となる、時刻（ｔ）を変数とする一次関数、または、この一次関数の平方根により増加係数を算出する。さらに、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が＋１より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０に増加係数を加算した値に変更してもよい。

また、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が－１より小さい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、－１以上、０に増加係数を加算した値以下における、Ｓ１０８で精鋭した乱数（ｒ）により定まる値に変更してもよい。さらに、更新部５０は、第１変数（ｘ_ｉ（ｔ＋Δｔ））が＋１より大きい場合、第１変数（ｘ_ｉ（ｔ＋Δｔ））を、０に増加係数を加算した値以上、＋１以下における、Ｓ１０８で生成した乱数（ｒ）により定まる値に変更してもよい。この場合、更新部５０は、式（２３）に含まれる指標値（ｘ_ａｖｅ）を増加係数に代えた式により、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を算出する。

続いて、Ｓ１２８において、更新部５０は、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））の制約処理をする。具体的には、更新部５０は、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））を、０、予め定められた値、または、乱数に応じた値に変更する。更新部５０は、乱数に応じた値に変更する場合、例えば、更新部５０は、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））を、予め定められた範囲（例えば、－０．１以上、＋０．１以下）内で一様な確率で発生する乱数に変更する。更新部５０は、Ｓ１１０を終えると、処理をＳ１１１に進める。

更新部５０は、以上のようなＳ１０５とＳ１１１との間のループ処理をＮ回実行することにより、次の処理を実行する。すなわち、更新部５０は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に基づき更新する。さらに、更新部５０は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が、第２値より大きい場合、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、符号を変更せずに、絶対値を０以上、第２値以下に変更する。さらに、更新部５０は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が、第２値より大きい場合、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））を、０、予め定められた値または乱数に応じた値に変更する。

更新部５０は、Ｓ１０５とＳ１１１との間のループ処理をＮ回実行した場合、処理をＳ１１２に進める。

続いて、更新部５０は、Ｓ１１２とＳ１１６との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ１１２とＳ１１６との間のループ処理において、更新部５０は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。

Ｓ１１３において、更新部５０は、Ｎ個の要素のそれぞれの対象時刻（ｔ＋Δｔ）における第１変数（ｘ_１（ｔ＋Δｔ）～ｘ_Ｎ（ｔ＋Δｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ＋Δｔ））を算出する。ＱＵＢＯ問題の場合、作用係数は、Ｊに含まれる結合係数およびｈに含まれる局所磁場係数である。ＨＯＢＯ問題の場合、作用係数は、Ｊ^（ｎ）に含まれる。

ＱＵＢＯ問題の場合、更新部５０は、式（２４）を算出する。

また、ＨＯＢＯ問題の場合、更新部５０は、式（２５）を算出する。

続いて、Ｓ１１４において、更新部５０は、対象時刻（ｔ＋Δｔ）における更新値（ｚ_ｉ（ｔ＋Δｔ））に係数（ｃ）を乗算することにより、外力（ｆ_ｉ（ｔ＋Δｔ））を算出する。具体的には、更新部５０は、式（２６）を算出する。

続いて、Ｓ１１５において、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に、外力（ｆ_ｉ（ｔ＋Δｔ））と対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））とに基づく値に単位時間（Δｔ）を乗算した値を、加算することにより、算出する。具体的には、更新部５０は、式（２７）を算出する。

更新部５０は、以上のようなＳ１１２とＳ１１６との間のループ処理をＮ回実行することにより、次の処理を実行する。すなわち、更新部５０は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、Ｎ個の要素のそれぞれの対象時刻（ｔ＋Δｔ）における第１変数（ｘ_１（ｔ＋Δｔ）～ｘ_Ｎ（ｔ＋Δｔ））に基づき更新する。

更新部５０は、Ｓ１１２とＳ１１６との間のループ処理をＮ回実行した場合、処理をＳ１１７に進める。Ｓ１１７において、更新部５０は、直前時刻（ｔ）に単位時間（Δｔ）を加算して、対象時刻（ｔ＋Δｔ）を更新する。Ｓ１１８において、更新部５０は、Ｓ１０４からＳ１１７までの処理を、ｔが終了時刻（Ｔ）を超えるまで繰り返す。そして、更新部５０は、ｔが終了時刻（Ｔ）より大きくなった場合、本フローを終了する。

そして、更新部５０は、Ｎ個の要素のそれぞれについて、終了時刻（ｔ＝Ｔ）における第１変数（ｘ_ｉ（Ｔ））の符号に応じて、対応するスピンの値を算出する。例えば、更新部５０は、終了時刻（ｔ＝Ｔ）における第１変数（ｘ_ｉ（Ｔ））の符号が負である場合、対応するスピンを－１とし、正である場合、対応するスピンを＋１とする。そして、更新部５０は、算出した複数のスピンの値を組合せ最適化問題の解として出力する。

図１８に示すフローチャートに従った処理を実行することにより、更新部５０は、Ｎ個の要素のそれぞれについて、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎の第１変数（ｘ_ｉ（ｔ＋Δｔ））および第２変数（ｙ_ｉ（ｔ＋Δｔ））を、単位時間毎に順次に、且つ、第１変数（ｘ_ｉ（ｔ＋Δｔ））および第２変数（ｙ_ｉ（ｔ＋Δｔ））を交互に更新する。また、図１８に示すフローチャートに従った処理を実行することにより、更新部５０は、単位時間毎に、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を算出した後に、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を算出する。これにより、更新部５０は、シンプレクティック・オイラー法を用いて第１アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

（第１アルゴリズムの処理フローの第２例）
図１９は、第１アルゴリズムを実行する場合における情報処理システム１００の処理の流れの第２例を示す。なお、図１９のフローチャートは、図１８で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。

情報処理システム１００は、式（９）および式（１０）に示す第１アルゴリズムを用いて最適化問題を解く場合、例えば、図１９に示す流れで処理を実行してもよい。

まず、Ｓ１０１およびＳ１０２において、更新部５０は、図１８に示す第１例と同様の処理を実行する。続いて、更新部５０は、Ｓ１０３とＳ１１８との間のループ処理を、図１８に示す第１例と同様に実行する。

Ｓ１０４において、更新部５０は、図１８に示す第１例と同様の処理を実行する。なお、更新部５０は、Ｓ１０４の処理をＳ１２６の直前で実行してもよい。

続いて、更新部５０は、Ｓ１２１とＳ１２５との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ１２１とＳ１２５との間のループ処理において、更新部５０は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。

Ｓ１２２において、更新部５０は、Ｎ個の要素のそれぞれの直前時刻（ｔ）における第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ））を算出する。

ＱＵＢＯ問題の場合、更新部５０は、式（２８）を算出する。

また、ＨＯＢＯ問題の場合、更新部５０は、式（２９）を算出する。

続いて、Ｓ１２３において、更新部５０は、直前時刻（ｔ）における更新値（ｚ_ｉ（ｔ））に係数（ｃ）を乗算することにより、直前時刻（ｔ）における外力（ｆ_ｉ（ｔ））を算出する。具体的には、更新部５０は、式（３０）を算出する。

続いて、Ｓ１２４において、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第２変数（ｙ_ｉ（ｔ））に、外力（ｆ_ｉ（ｔ））と対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））とに基づく値に単位時間（Δｔ）を乗算した値を、加算することにより、算出する。具体的には、更新部５０は、式（３１）を算出する。

更新部５０は、以上のようなＳ１２１とＳ１２５との間のループ処理をＮ回実行することにより、次の処理を実行する。すなわち、更新部５０は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、Ｎ個の要素のそれぞれの直前時刻（ｔ）における第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））に基づき更新する。

更新部５０は、Ｓ１２１とＳ１２５との間のループ処理をＮ回実行した場合、処理をＳ１２６に進める。

続いて、更新部５０は、Ｓ１２６とＳ１２９との間のループ処理を、ｉ＝１からｉ＝Ｎまでｉを１ずつインクリメントしながら繰り返す。Ｓ１２６とＳ１２９との間のループ処理において、更新部５０は、Ｎ個の要素のうちのｉ番目の要素を、対象要素として処理を実行する。

Ｓ１２７において、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の直前時刻（ｔ）における第１変数（ｘ_ｉ（ｔ））に、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））と予め定められた定数（Ｄ）と単位時間（Δｔ）とを乗算した値を加算することにより算出する。具体的には、更新部５０は、式（３２）を算出する。

続いて、Ｓ１０７において、更新部５０は、図１８に示す第１例と同様の処理を実行する。更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が第２値以下である場合、処理をＳ１２９に進める（Ｓ１０７のＮｏ）。更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が第２値より大きい場合、処理をＳ１０８に進める。

続いて、Ｓ１０８およびＳ１０９において、更新部５０は、図１８に示す第１例と同様の処理を実行する。

続いて、Ｓ１１０において、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））の制約処理をする。具体的には、更新部５０は、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、０、予め定められた値、または、乱数に応じた値に変更する。更新部５０は、乱数に応じた値に変更する場合、例えば、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、予め定められた範囲（例えば、－０．１以上、＋０．１以下）内で一様な確率で発生する乱数に変更する。更新部５０は、Ｓ１２８を終えると、処理をＳ１２９に進める。

更新部５０は、以上のようなＳ１２６とＳ１２９との間のループ処理をＮ回実行することにより、次の処理を実行する。すなわち、更新部５０は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））に基づき更新する。さらに、更新部５０は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が、第２値より大きい場合、対象要素の対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を、符号を変更せずに、絶対値を０以上、第２値以下に変更する。さらに、更新部５０は、Ｎ個の要素のそれぞれについて、対象要素の対象時刻（ｔ＋Δｔ）における第１変数の絶対値（｜ｘ_ｉ（ｔ＋Δｔ）｜）が、第２値より大きい場合、対象要素の対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を、０、予め定められた値または乱数に応じた値に変更する。

更新部５０は、Ｓ１２６とＳ１２９との間のループ処理をＮ回実行した場合、処理をＳ１１７に進める。Ｓ１１７において、更新部５０は、図１８に示す第１例と同様の処理を実行する。Ｓ１１８において、更新部５０は、Ｓ１０４からＳ１１７までの処理を、ｔが終了時刻（Ｔ）を超えるまで繰り返す。そして、更新部５０は、ｔが終了時刻（Ｔ）より大きくなった場合、本フローを終了する。さらに、更新部５０は、Ｎ個の要素のそれぞれについて、終了時刻（ｔ＝Ｔ）における第１変数（ｘ_ｉ（Ｔ））の符号に応じて、対応するスピンの値を算出し、算出した複数のスピンの値を組合せ最適化問題の解として出力する。

図１９に示すフローチャートに従った処理を実行することにより、更新部５０は、Ｎ個の要素のそれぞれについて、初期時刻（ｔ＝０）から終了時刻（ｔ＝Ｔ）まで単位時間（Δｔ）毎の第１変数（ｘ_ｉ（ｔ＋Δｔ））および第２変数（ｙ_ｉ（ｔ＋Δｔ））を、単位時間毎に順次に、且つ、第１変数（ｘ_ｉ（ｔ＋Δｔ））および第２変数（ｙ_ｉ（ｔ＋Δｔ））を交互に更新する。また、図１９に示すフローチャートに従った処理を実行することにより、更新部５０は、単位時間毎に、対象時刻（ｔ＋Δｔ）における第２変数（ｙ_ｉ（ｔ＋Δｔ））を算出した後に、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｉ（ｔ＋Δｔ））を算出する。これにより、更新部５０は、シンプレクティック・オイラー法を用いて第１アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

（第２アルゴリズムの処理フローの第１例）
図２０は、第２アルゴリズムを実行する場合における情報処理システム１００の処理の流れの第１例を示す。

情報処理システム１００は、式（１１）および式（１２）に示す第２アルゴリズムを用いて最適化問題を解く場合、例えば、図２０に示す流れで処理を実行する。なお、図２０のフローチャートは、図１８で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図１８で示した流れで第１アルゴリズムを実行する場合との相違点を説明する。

図２０の処理は、Ｓ１１１とＳ１１２との間において、Ｓ２０１の処理が追加されている。Ｓ２０１において、更新部５０は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ＋Δｔ）における第１変数（ｘ_ｊ（ｔ＋Δｔ））の符号（ｓ_ｊ（ｔ＋Δｔ））を算出する。具体的には、ｊ＝１～Ｊ＝Ｎまでのそれぞれについて、更新部５０は、式（３３）を算出する。

また、更新部５０は、Ｓ１１３に代えて、Ｓ２０２を実行する。Ｓ２０２において、更新部５０は、Ｎ個の要素のそれぞれの対象時刻（ｔ＋Δｔ）における第１変数の符号（ｓ_１（ｔ＋Δｔ）～ｓ_Ｎ（ｔ＋Δｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ＋Δｔ））を算出する。

ＱＵＢＯ問題の場合、更新部５０は、式（３４）を算出する。

また、ＨＯＢＯ問題の場合、更新部５０は、式（３５）を算出する。

図２０に示すフローチャートに従った処理を実行することにより、更新部５０は、シンプレクティック・オイラー法を用いて第２アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

（第２アルゴリズムの処理フローの第２例）
図２１は、第２アルゴリズムを実行する場合における情報処理システム１００の処理の流れの第２例を示す。

情報処理システム１００は、式（１１）および式（１２）に示す第２アルゴリズムを用いて最適化問題を解く場合、例えば、図２１に示す流れで処理を実行することもできる。なお、図２１のフローチャートは、図１９で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図１９で示した流れで第１アルゴリズムを実行する場合との相違点を説明する。

図２１の処理は、Ｓ１２１の前において、Ｓ２１１の処理が追加されている。なお、Ｓ２１１は、Ｓ１０４の前に追加されてもよい。

Ｓ２１１において、更新部５０は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ）における第１変数（ｘ_ｊ（ｔ））の符号（ｓ_ｊ（ｔ））を算出する。更新部５０は、式（３６）を演算することにより、Ｎ個の要素のそれぞれについて、直前時刻（ｔ）における第１変数（ｘ_ｊ（ｔ））の符号（ｓ_ｊ（ｔ））を算出する。具体的には、ｊ＝１～Ｊ＝Ｎまでのそれぞれについて、更新部５０は、式（３６）を算出する。

また、更新部５０は、Ｓ１２２に代えて、Ｓ２１２を実行する。Ｓ２１２において、更新部５０は、Ｎ個の要素のそれぞれの直前時刻（ｔ）における第１変数の符号（ｓ_１（ｔ）～ｓ_Ｎ（ｔ））と、対象要素とＮ個の要素のそれぞれとの組毎に予め定められる作用係数と、に基づき更新値（ｚ_ｉ（ｔ））を算出する。

ＱＵＢＯ問題の場合、更新部５０は、式（３７）を算出する。

また、ＨＯＢＯ問題の場合、更新部５０は、式（３８）を算出する。

図２１に示すフローチャートに従った処理を実行することにより、更新部５０は、シンプレクティック・オイラー法を用いて第２アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

（第３アルゴリズムの処理フローの第１例）
図２２は、第３アルゴリズムを実行する場合における情報処理システム１００の処理の流れの第１例を示す。

情報処理システム１００は、式（１３）および式（１４）に示す第３アルゴリズムを用いて最適化問題を解く場合、例えば、図２２に示す流れで処理を実行する。なお、図２２のフローチャートは、図１８で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図１８で示した流れで第１アルゴリズムを実行する場合との相違点を説明する。

更新部５０は、Ｓ１１４に代えて、Ｓ３０１を実行する。Ｓ３０１において、更新部５０は、対象時刻（ｔ＋Δｔ）における更新値（ｚ_ｉ（ｔ＋Δｔ））の符号に、予め定められた関数により定まる係数を乗算することにより、対象時刻（ｔ＋Δｔ）における外力（ｆ_ｉ（ｔ＋Δｔ））を算出する。具体的には、更新部５０は、式（３９）を算出する。

なお、ｇ（ｔ＋Δｔ）は、式（４０）により表される。

図２２に示すフローチャートに従った処理を実行することにより、更新部５０は、シンプレクティック・オイラー法を用いて第３アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

（第３アルゴリズムの処理フローの第２例）
図２３は、第３アルゴリズムを実行する場合における情報処理システム１００の処理の流れの第２例を示す。

情報処理システム１００は、式（１３）および式（１４）に示す第３アルゴリズムを用いて最適化問題を解く場合、例えば、図２３に示す流れで処理を実行することもできる。なお、図２３のフローチャートは、図１９で示したフローチャートの処理と同一の処理のステップには、同一のステップ番号を付けている。以下、図１９で示した流れで第１アルゴリズムを実行する場合との相違点を説明する。

更新部５０は、Ｓ１２３に代えて、Ｓ３１１を実行する。Ｓ３１１において、更新部５０は、直前時刻（ｔ）における更新値（ｚ_ｉ（ｔ））の符号に、予め定められた関数により定まる係数を乗算することにより、直前時刻（ｔ）における外力（ｆ_ｉ（ｔ））を算出する。具体的には、更新部５０は、式（４１）を演算することにより、直前時刻（ｔ）における外力（ｆ_ｉ（ｔ））を算出する。

なお、ｇ（ｔ）は、式（４２）により表される。

図２３に示すフローチャートに従った処理を実行することにより、更新部５０は、シンプレクティック・オイラー法を用いて第３アルゴリズムに従った演算を実行して、終了時刻（ｔ＝Ｔ）におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））およびＮ個の第２変数（ｙ_１（ｔ）～ｙ_Ｎ（ｔ））を算出することができる。

（半導体装置への回路化）
図２４は、式（９）の第１アルゴリズム、式（１１）の第２アルゴリズムまたは式（１３）の第３アルゴリズムを実行することにより、組合せ最適化問題の解を算出する演算装置６０の構成を示す図である。演算装置６０は、ＦＰＧＡ、ゲートアレイまたは特定用途向け集積回路等であって、半導体装置に回路化して実現される。

演算装置６０は、演算回路６１と、入力回路６２と、出力回路６３と、設定回路６４とを備える。

演算回路６１は、時刻を表すパラメータであるｔを、初期時刻（例えば０）から終了時刻まで、単位時間（Δｔ）ずつ増加させる。演算回路６１は、Ｎ個の要素（仮想的な粒子）のそれぞれに対応付けられた、第１変数（ｘ_ｉ）および第２変数（ｙ_ｉ）を算出する。第１変数（ｘ_ｉ）は、対応する要素（仮想的な粒子）の位置を表す。第２変数（ｙ_ｉ）は、対応する要素（仮想的な粒子）の運動量を表す。

演算回路６１は、初期時刻から終了時刻まで単位時間毎に、Ｎ個の第１変数（ｘ_ｉ）のそれぞれおよびＮ個の第２変数（ｙ_ｉ）のそれぞれを、単位時間毎に順次に且つ第１変数（ｘ_ｉ）と第２変数（ｙ_ｉ）とを交互に算出する。より具体的には、演算回路６１は、初期時刻から終了時刻まで単位時間毎に、式（９）、式（１１）または式（１３）のアルゴリズムで示された演算を実行する。そして、演算回路６１は、終了時刻におけるＮ個の第１変数（ｘ_ｉ）のそれぞれの値（すなわち、Ｎ個の仮想的な粒子のそれぞれの位置）を２値化することにより、組合せ最適化問題の解を算出する。

入力回路６２は、演算回路６１による演算処理に先だって、初期時刻におけるＮ個の第１変数（ｘ_ｉ）およびＮ個の第２変数（ｙ_ｉ）のそれぞれの初期値（すなわち、複数の仮想的な粒子のそれぞれの初期位置および初期運動量）を取得して、演算回路６１に与える。出力回路６３は、演算回路６１による演算処理の終了後に、組合せ最適化問題の解を演算回路６１から取得する。そして、出力回路６３は、取得した解を出力する。設定回路６４は、演算回路６１による演算処理に先だって、各パラメータを演算回路６１に対して設定する。

（第１アルゴリズムの処理を実行する演算回路６１）
図２５は、式（９）の第１アルゴリズムを実行する演算回路６１のブロック構成を示す図である。演算回路６１は、式（９）の第１アルゴリズムを実行する場合、図２５に示すように構成される。なお、図２５の説明において、直前時刻は、ｔ_１と表される。また、対象時刻は、ｔ_２と表される。図２６および図２７でも同様である。

式（９）の第１アルゴリズムを実行する演算回路６１は、Ｘメモリ６６と、Ｙメモリ６７と、作用演算回路６８と、更新回路６９と、制御回路７０とを備える。

Ｘメモリ６６は、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｉ（ｔ_１））を記憶する。Ｘメモリ６６は、時刻の更新に伴って、記憶している直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｉ（ｔ_１））が書き換えられる。すなわち、Ｘメモリ６６は、対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｉ（ｔ_２））が算出された場合、算出された対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｉ（ｔ_２））が、新たな直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｉ（ｔ_１））として書き込まれる。設定回路６４は、演算に先だって、初期時刻におけるＮ個の第１変数ｘ_ｉをＸメモリ６６に書き込む。

Ｙメモリ６７は、直前時刻（ｔ_１）におけるＮ個の第２変数（ｙ_ｉ（ｔ_１））を記憶する。Ｙメモリ６７は、時刻の更新に伴って、記憶している直前時刻（ｔ_１）におけるＮ個の第２変数（ｙ_ｉ（ｔ_１））が書き換えられる。すなわち、Ｙメモリ６７は、対象時刻（ｔ_２）におけるＮ個の第２変数（ｙ_ｉ（ｔ_２））が算出された場合、算出された対象時刻（ｔ_２）におけるＮ個の第２変数（ｙ_ｉ（ｔ_２））が、新たな直前時刻（ｔ_１）におけるＮ個の第２変数（ｙ_ｉ（ｔ_１））として書き込まれる。設定回路６４は、演算に先だって、初期時刻におけるＮ個の第２変数ｙ_ｉをＹメモリ６７に書き込む。

作用演算回路６８は、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））をＸメモリ６６から取得する。そして、作用演算回路６８は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における更新値（ｚ_ｉ（ｔ_１））を算出する。

更新回路６９は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における更新値（ｚ_ｉ（ｔ_１））を作用演算回路６８から取得する。さらに、更新回路６９は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））をＸメモリ６６から取得し、直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））をＹメモリ６７から取得する。そして、更新回路６９は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出し、Ｘメモリ６６に記憶されている直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））を書き換える。これとともに、更新回路６９は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出し、Ｙメモリ６７に記憶されている直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））を書き換える。

制御回路７０は、対象時刻（ｔ_２）を単位時間（Δｔ）毎に順次に更新することにより、作用演算回路６８および更新回路６９に対して単位時間（Δｔ）毎の第１変数（ｘ_ｉ（ｔ））および第２変数（ｙ_ｉ（ｔ））を順次に算出させる。

さらに、制御回路７０は、１からＮまでのインクリメントすることによりインデックス（ｉ）を発生して、作用演算回路６８および更新回路６９に対して、Ｎ個の要素のそれぞれに対応する、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））および対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））をインデックス順に算出させる。なお、作用演算回路６８および更新回路６９は、複数個のインデックスに対応する複数個の第１変数（ｘ_ｉ（ｔ_２））および複数個の第２変数（ｙ_ｉ（ｔ_２））を並列に算出してもよい。

作用演算回路６８は、Ｊメモリ７１と、Ｈメモリ７２と、行列演算回路７３と、α関数回路７４と、第１加算回路７５とを有する。

Ｊメモリ７１は、（Ｎ×Ｎ）個の結合係数を含むＮ×Ｎの行列を記憶する。Ｊ_ｉ，_ｊは、行列に含まれるｉ行ｊ列の結合係数を表す。Ｊ_ｉ，_ｊは、組合せ最適化問題を表すイジングモデルにおける、ｉ番目のスピンとｊ番目のスピンとの結合係数を表す。設定回路６４は、演算に先だって、ユーザ等により予め生成された行列をＪメモリ７１に書き込む。

Ｈメモリ７２は、Ｎ個の局所磁場係数を含む配列を記憶する。ｈ_ｉは、配列に含まれるｉ番目の局所磁場係数を表す。ｈ_ｉは、組合せ最適化問題を表すイジングモデルにおける、ｉ番目のスピンに作用する局所磁場を表す。設定回路６４は、演算に先だって、ユーザ等により予め生成された配列をＨメモリ７２に書き込む。

行列演算回路７３は、Ｘメモリ６６から、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））を取得する。行列演算回路７３は、Ｎ個の要素のそれぞれについて、Ｊメモリ７１から対象の行に含まれるＮ個の結合係数Ｊ_ｉ，_ｊを取得する。そして、行列演算回路７３は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））と、対象の行に含まれるＮ個の結合係数Ｊ_ｉ，_ｊとの積和演算を実行する。

α関数回路７４は、Ｎ個の要素のそれぞれについて、Ｈメモリ７２から対象の局所磁場係数ｈ_ｉを取得する。α関数回路７４は、Ｎ個の要素のそれぞれについて、｛－ｈ_ｉα（ｔ_１）｝の演算を実行する。α（ｔ）は、予め設定された関数である。

第１加算回路７５は、Ｎ個の要素のそれぞれについて、行列演算回路７３による積和演算結果と、α関数回路７４による演算結果とを加算する。この演算により、第１加算回路７５は、Ｎ個の要素のそれぞれについて、式（４３）で表される直前時刻（ｔ_１）における更新値（ｚ_ｉ（ｔ_１））を出力する。

更新回路６９は、第１乗算回路７９と、Ｐ関数回路８０と、第２乗算回路８１と、第２加算回路８２と、第３乗算回路８３と、第３加算回路８４と、制約前Ｙメモリ８５と、第４乗算回路８６と、第４加算回路８７と、平均化回路９０と、判定回路９１と、Ｘ制約回路９２と、Ｙ制約回路９３とを有する。

第１乗算回路７９は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における更新値（ｚ_ｉ（ｔ_１））に係数である－ｃを乗算する。Ｐ関数回路８０は、Ｎ個の要素のそれぞれについて、｛－Ｄ＋ｐ（ｔ_１）｝の演算を実行する。第２乗算回路８１は、Ｎ個の要素のそれぞれについて、Ｘメモリ６６から直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））を取得する。そして、第２乗算回路８１は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））と、Ｐ関数回路８０の演算結果とを乗算する。

第２加算回路８２は、Ｎ個の要素のそれぞれについて、第１乗算回路７９の演算結果と、第２乗算回路８１の演算結果とを加算する。第３乗算回路８３は、第２加算回路８２の演算結果に単位時間であるΔｔを乗算する。

第３加算回路８４は、Ｎ個の要素のそれぞれについて、Ｙメモリ６７から直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））を取得する。第３加算回路８４は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第２変数（ｙ_ｉ（ｔ_１））と、第３乗算回路８３の演算結果とを加算する。この演算により、第３加算回路８４は、Ｎ個の要素のそれぞれについて、式（４４）で表される対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を出力する。

そして、第３加算回路８４は、算出したＮ個の要素のそれぞれの対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を、制約前Ｙメモリ８５に書き込む。制約前Ｙメモリ８５は、Ｙ制約回路９３による制約前の、対象時刻（ｔ_２）におけるＮ個の第２変数（ｙ_ｉ（ｔ_２））を記憶する。

第４乗算回路８６は、Ｎ個の要素のそれぞれについて、制約前Ｙメモリ８５から対象時刻（ｔ_２）における、制約前の第２変数（ｙ_ｉ（ｔ_２））を取得する。第４乗算回路８６は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における制約前の第２変数（ｙ_ｉ（ｔ_２））に、｛ＤΔｔ｝を乗算する。

第４加算回路８７は、Ｎ個の要素のそれぞれについて、Ｘメモリ６６から直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））を取得する。第４加算回路８７は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における第１変数（ｘ_ｉ（ｔ_１））と、第４乗算回路８６の演算結果とを加算する。この演算により、第４加算回路８７は、Ｎ個の要素のそれぞれについて、式（４５）で表される対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を出力する。

平均化回路９０は、Ｘメモリ６６に記憶された直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_１（ｔ_１）～ｘ_Ｎ（ｔ_１））の大きさの平均を表す指標値（ｘ_ａｖｅ）を算出する。例えば、指標値（ｘ_ａｖｅ）は、直前時刻におけるＮ個の第１変数（ｘ_１（ｔ）～ｘ_Ｎ（ｔ））の二乗平均平方根または絶対値平均である。

判定回路９１は、Ｎ個の要素のそれぞれについて、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数の絶対値（｜ｘ_ｉ（ｔ_２）｜）が、予め定められた第２値より大きいか否かを判断する。例えば、第２値は、＋１である。判定回路９１は、対象時刻（ｔ_２）における第１変数の絶対値（｜ｘ_ｉ（ｔ_２）｜）が第２値より大きい場合、Ｘ制約回路９２およびＹ制約回路９３にイネーブル信号（ＥＮ）を与える。

Ｘ制約回路９２は、Ｎ個の要素のそれぞれについて、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を受け取る。Ｘ制約回路９２は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取らなかった場合、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を、そのままＸメモリ６６に書き込む。

Ｘ制約回路９２は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取った場合、第４加算回路８７により算出された対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））に対して制約処理を実行し、制約処理をした後の第１変数（ｘ_ｉ（ｔ_２））をＸメモリ６６に書き込む。

ここで、Ｘ制約回路９２は、Ｎ個の要素のそれぞれについて、制約処理として、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を、符号を変更せずに、絶対値を０以上、第２値（例えば＋１）以下の値に変更する。

例えば、Ｎ個の要素のそれぞれについて、Ｘ制約回路９２は、第１変数の絶対値（｜ｘ_ｉ（ｔ_２）｜）を、予め定められた値、または、０以上第２値（例えば＋１）以下の範囲における所定区間で一様な確率で発生する乱数に応じた値に設定する。この場合、所定区間は、０以上、第２値以下の範囲の内側であれば、どのような範囲であってもよい。

また、Ｎ個の要素のそれぞれについて、Ｘ制約回路９２は、象時刻（ｔ_２）における第１変数の絶対値（｜ｘ_ｉ（ｔ_２）｜）を、平均化回路９０で算出した指標値に変更してもよい。

また、Ｎ個の要素のそれぞれについて、Ｘ制約回路９２は、対象時刻（ｔ_２）における第１変数の絶対値（｜ｘ_ｉ（ｔ_２）｜）を、平均化回路９０で算出した指標値以上、第２値（例えば＋１）以下において、乱数により定まる値に変更してもよい。この場合、乱数により定まる値は、例えば、指標値と第２値（例えば＋１）以下の範囲内の所定区間で一様な確率で発生する乱数である。

また、Ｎ個の要素のそれぞれについて、Ｘ制約回路９２は、対象時刻（ｔ_２）における第１変数の絶対値（｜ｘ_ｉ（ｔ_２）｜）を、０から第２値以下の間で初期時刻から終了時刻まで時間経過に従って増加する増加係数に変更してもよい。この場合、増加関数は、例えば、初期時刻で０、終了時刻で第２値となる時刻を変数とする一次関数、または、この一次関数の平方根である。また、Ｎ個の要素のそれぞれについて、Ｘ制約回路９２は、対象時刻（ｔ_２）における第１変数の絶対値（｜ｘ_ｉ（ｔ_２）｜）を、増加係数以上、第２値（例えば＋１）以下における、乱数により定まる値に変更してもよい。この場合、乱数により定まる値は、例えば、増加係数と第２値（例えば＋１）との間で一様な確率で発生する乱数である。

Ｙ制約回路９３は、Ｎ個の要素のそれぞれについて、制約前Ｙメモリ８５から、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を取得する。Ｙ制約回路９３は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取らなかった場合、制約前Ｙメモリ８５から取得した対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を、そのままＹメモリ６７に書き込む。

Ｙ制約回路９３は、Ｎ個の要素のそれぞれについて、判定回路９１からイネーブル信号（ＥＮ）を受け取った場合、第４加算回路８７により算出された対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））に対して制約処理を実行し、制約処理をした後の第２変数（ｙ_ｉ（ｔ_２））をＹメモリ６７に書き込む。

ここで、Ｙ制約回路９３は、Ｎ個の要素のそれぞれについて、制約処理として、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を、第２変数（ｙ_ｉ（ｔ_２））に乱数を乗じた値、０、予め定められた値、または、乱数に応じた値に変更する。なお、この場合において、Ｙ制約回路９３は、変更後の値が、所定の範囲内の値とするように処理を行う。例えば、Ｙ制約回路９３は、変更後の値が、－０．１以上、＋０．１以下の範囲内となるように処理を行ってもよい。

以上により、演算回路６１は、式（９）の第１アルゴリズムを実行して、終了時刻（Ｔ）におけるＮ個の第１変数（ｘ_ｉ（Ｔ））およびＮ個の第２変数（ｙ_ｉ（Ｔ））を算出することができる。さらに、演算回路６１は、更新によって、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_１））の絶対値が第２値より大きくなった場合、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_１））を、符号を保持したまま、絶対値を０以上第２値以下の値に変更することができる。

なお、上述した演算回路６１は、単位時間毎に、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出した後に、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出する。これに代えて、演算回路６１は、単位時間毎に、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出した後に、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出してもよい。

この場合、制約前Ｙメモリ８５は、Ｙ制約回路９３による制約前の、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｉ（ｔ_１））を記憶する。そして、第４加算回路８７は、Ｎ個の要素のそれぞれについて、式（４６）で表される対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を出力する。

行列演算回路７３は、Ｘメモリ６６から、対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｊ（ｔ_２））を取得する。そして、第１加算回路７５は、Ｎ個の要素のそれぞれについて、式（４７）で表される対象時刻（ｔ_２）における更新値（ｚ_ｉ（ｔ_２））を出力する。

また、第１乗算回路７９は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における更新値（ｚ_ｉ（ｔ_２））に－ｃを乗算する。Ｐ関数回路８０は、Ｎ個の要素のそれぞれについて、｛－Ｄ＋ｐ（ｔ_２）｝の演算を実行する。第２乗算回路８１は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））と、Ｐ関数回路８０の演算結果とを乗算する。

そして、第３加算回路８４は、Ｎ個の要素のそれぞれについて、式（４８）で表される対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を出力する。

このような処理によっても、演算回路６１は、式（９）の第１アルゴリズムを実行して、終了時刻（Ｔ）におけるＮ個の第１変数（ｘ_ｉ（Ｔ））およびＮ個の第２変数（ｙ_ｉ（Ｔ））を算出することができる。さらに、演算回路６１は、更新によって、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））の絶対値が第２値より大きくなった場合、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を、符号を保持したまま、絶対値を０以上第２値以下の値に変更することができる。

（第２アルゴリズムの処理を実行する演算回路６１）
図２６は、式（１１）の第２アルゴリズムを実行する演算回路６１のブロック構成を示す図である。演算回路６１は、式（１１）の第２アルゴリズムを実行する場合、図２５に示すように構成される。なお、図２６に示す演算回路６１は、図２５に示した構成とほぼ同一の構成を有するので、同一の機能の構成要素については同一の符号を付けて詳細な説明を省略する。

式（１１）の第２アルゴリズムを実行する演算回路６１は、図２５に示した構成と比較して、符号化回路９６をさらに有する。

符号化回路９６は、Ｘメモリ６６から直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））のそれぞれを取得する。符号化回路９６は、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））のそれぞれの符号（－１または＋１）を抽出し、Ｎ個の第１変数の符号（ｓ_ｊ（ｔ_１））を出力する。

行列演算回路７３は、直前時刻（ｔ_１）におけるＮ個の第１変数（ｘ_ｊ（ｔ_１））に代えて、符号化回路９６から、Ｎ個の第１変数の符号（ｓ_ｊ（ｔ_１））を取得する。そして、行列演算回路７３は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）におけるＮ個の第１変数の符号（ｓ_ｊ（ｔ_１））と、対象の行に含まれるＮ個の結合係数Ｊ_ｉ，_ｊとの積和演算を実行する。

以上により、演算回路６１は、式（１１）の第２アルゴリズムを実行して、終了時刻（Ｔ）におけるＮ個の第１変数（ｘ_ｉ（Ｔ））およびＮ個の第２変数（ｙ_ｉ（Ｔ））を算出することができる。

なお、演算回路６１が、単位時間毎に、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出した後に、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出する場合、符号化回路９６は、Ｘメモリ６６から対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｊ（ｔ_２））のそれぞれを取得する。符号化回路９６は、対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｊ（ｔ_２））のそれぞれの符号（－１または＋１）を抽出し、Ｎ個の第１変数の符号（ｓ_ｊ（ｔ_２））を出力する。この場合、行列演算回路７３は、対象時刻（ｔ_２）におけるＮ個の第１変数（ｘ_ｊ（ｔ_２））に代えて、符号化回路９６から、Ｎ個の第１変数の符号（ｓ_ｊ（ｔ_２））を取得する。そして、行列演算回路７３は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）におけるＮ個の第１変数の符号（ｓ_ｊ（ｔ_２））と、対象の行に含まれるＮ個の結合係数Ｊ_ｉ，_ｊとの積和演算を実行する。このような処理をしても、演算回路６１は、式（１１）の第２アルゴリズムを実行して、終了時刻（Ｔ）におけるＮ個の第１変数（ｘ_ｉ（Ｔ））およびＮ個の第２変数（ｙ_ｉ（Ｔ））を算出することができる。

（第３アルゴリズムの処理を実行する演算回路６１）
図２７は、式（１３）の第３アルゴリズムを実行する演算回路６１のブロック構成を示す図である。演算回路６１は、式（１３）の第３アルゴリズムを実行する場合、図２７に示すように構成される。なお、図２７に示す演算回路６１は、図２５に示した構成とほぼ同一の構成を有するので、同一の機能の構成要素については同一の符号を付けて詳細な説明を省略する。

式（１３）の第３アルゴリズムを実行する演算回路６１は、図２５に示した構成と比較して、更新回路６９の内部構成が異なる。式（１３）の第３アルゴリズムを実行する更新回路６９は、図２５に示した構成と比較して、符号化回路９６と、Ｇ関数回路９７とをさらに有する。

符号化回路９６は、Ｎ個の要素のそれぞれについて、作用演算回路６８から、直前時刻（ｔ_１）における更新値（ｚ_ｉ（ｔ_１））を取得する。符号化回路９６は、Ｎ個の要素のそれぞれについて、直前時刻（ｔ_１）における更新値（ｚ_ｉ（ｔ_１））から符号（－１または＋１）を抽出し、直前時刻（ｔ_１）における更新値の符号｛ｓｇｎ（ｚ_ｉ（ｔ_１））｝を出力する。

Ｇ関数回路９７は、Ｎ個の要素のそれぞれについて、予め定められた関数ｇ（ｔ）の演算を実行する。具体的には、Ｇ関数回路９７は、ｇ（ｔ_１）＝｛Ｄ－ｐ（ｔ_１）｝√（ｐ（ｔ_１））の演算を実行する。

そして、第１乗算回路７９は、Ｎ個の要素のそれぞれについて、符号化回路９６から出力された直前時刻（ｔ_１）における更新値の符号｛ｓｇｎ（ｚ_ｉ（ｔ_１））｝と、Ｇ関数回路９７の演算結果とを乗算する。

以上により、演算回路６１は、式（１３）の第３アルゴリズムを実行して、終了時刻（Ｔ）におけるＮ個の第１変数（ｘ_ｉ（Ｔ））およびＮ個の第２変数（ｙ_ｉ（Ｔ））を算出することができる。

なお、演算回路６１が、単位時間毎に、対象時刻（ｔ_２）における第１変数（ｘ_ｉ（ｔ_２））を算出した後に、対象時刻（ｔ_２）における第２変数（ｙ_ｉ（ｔ_２））を算出する場合、符号化回路９６は、作用演算回路６８から、対象時刻（ｔ_２）における更新値（ｚ_ｉ（ｔ_２））を取得する。符号化回路９６は、Ｎ個の要素のそれぞれについて、対象時刻（ｔ_２）における更新値（ｚ_ｉ（ｔ_２））から符号（－１または＋１）を抽出し、対象時刻（ｔ_２）における更新値の符号｛ｓｇｎ（ｚ_ｉ（ｔ_２））｝を出力する。この場合、Ｇ関数回路９７は、Ｎ個の要素のそれぞれについて、ｇ（ｔ_２）＝｛Ｄ－ｐ（ｔ_２）｝√（ｐ（ｔ_２））の演算を実行する。そして、第１乗算回路７９は、Ｎ個の要素のそれぞれについて、符号化回路９６から出力された対象時刻（ｔ_２）における更新値の符号｛ｓｇｎ（ｚ_ｉ（ｔ_２））｝と、Ｇ関数回路９７の演算結果とを乗算する。

なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。

１ 管理サーバ
２ ネットワーク
３、３ａ、３ｂ、３ｃ 計算サーバ
４、４ａ、４ｂ、４ｃ ケーブル
５ スイッチ
６ 情報端末
１０ プロセッサ
１１ 管理部
１２ 変換部
１３ 制御部
１４ 記憶部
１４Ａ 問題データ
１４Ｂ 計算データ
１４Ｃ 管理プログラム
１４Ｄ 変換プログラム
１４Ｅ、３４Ｃ 制御プログラム
１５、３１ 通信回路
１６ 入力回路
１７ 出力回路
１８ 操作装置
１９ 表示装置
２０ バス
３２ 共有メモリ
３３Ａ、３３Ｂ、３３Ｃ、３３Ｄ プロセッサ
３４ ストレージ
３４Ａ 計算データ
３４Ｂ 計算プログラム
３５ ホストバスアダプタ
５１ 作用演算部
１００ 情報処理システム

Claims (22)

1. 組合せ最適化問題を解く計算装置であって、
   第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎の前記第１変数および前記第２変数を、前記単位時間毎に順次に且つ前記第１変数および前記第２変数を交互に更新する更新部と、
   前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力部と、
   を備え、
   前記複数の要素は、前記組合せ最適化問題を表す複数の離散変数に対応し、
   前記第１変数および前記第２変数のそれぞれは、実数により表され、
   前記単位時間毎の算出処理において、前記更新部は、前記複数の要素のそれぞれについて、
   前記第１変数を前記第２変数に基づき更新し、
   前記第２変数を前記複数の要素のそれぞれの前記第１変数に基づき更新し、
   前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値以上、予め定められたしきい値以下の値に変更し、
   前記第１変数が予め定められた第２値より大きい場合、前記第１変数を、前記しきい値以上、前記第２値以下の値に変更し、
   前記第２値は、前記第１値より大きく、
   前記しきい値は、前記第１値より大きく前記第２値より小さい
   計算装置。
2. 前記出力部は、
   前記終了時刻における前記複数の要素のそれぞれについて、前記第１変数を前記しきい値により２値化した離散変数の値を算出し、
   算出した前記複数の離散変数の値を前記組合せ最適化問題の解として出力する
   請求項１に記載の計算装置。
3. 前記単位時間毎の算出処理において、前記更新部は、
   前記複数の要素のそれぞれについて、
   前記第１変数が前記第１値より小さい場合、前記第１変数を、前記第１値以上前記しきい値以下の予め定められた値、または、前記第１値以上前記しきい値以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値に変更し、
   前記第１変数が前記第２値より大きい場合、前記第１変数を、前記しきい値以上前記第２値以下の予め定められた値、または、前記しきい値以上前記第２値以下の範囲内における所定区間で一様な確率で発生する乱数に応じた値に変更する
   請求項１または２に記載の計算装置。
4. 前記単位時間毎の算出処理において、前記更新部は、
   前記複数の要素のそれぞれについて、
   前記第１変数が前記第１値より小さい場合、前記第１変数を、前記しきい値から指標値を減算した値に変更し、
   前記第１変数が前記第２値より大きい場合、前記第１変数を、前記しきい値に前記指標値を加算した値に変更し、
   前記指標値は、前記複数の要素の前記第１変数の前記しきい値からの偏差の大きさに対する平均を表す
   請求項１または２に記載の計算装置。
5. 前記単位時間毎の算出処理において、前記更新部は、
   前記複数の要素のそれぞれについて、
   前記第１変数が前記第１値より小さい場合、前記第１変数を、前記第１値以上、前記しきい値から指標値を減じた値以上の、乱数により定まる値に変更し、
   前記第１変数が前記第２値より大きい場合、前記第１変数を、前記しきい値に前記指標値を加算した値以上、前記第２値以下の、乱数により定まる値に変更し、
   前記指標値は、前記複数の要素の前記第１変数の前記しきい値からの偏差の大きさに対する平均を表す
   請求項１または２に記載の計算装置。
6. 前記指標値は、前記複数の要素の前記第１変数の前記しきい値からの偏差の二乗平均平方根または絶対値平均である
   請求項４または５に記載の計算装置。
7. 前記単位時間毎の算出処理において、前記更新部は、
   前記複数の要素のそれぞれについて、
   前記第１変数が前記第１値より小さい場合、前記第１変数を、前記しきい値から増加係数を減算した値に変更し、
   前記第１変数が前記第２値より大きい場合、前記第１変数を、前記しきい値に前記増加係数を加算した値に変更し、
   前記増加係数は、前記初期時刻において０であり、前記初期時刻から前記終了時刻まで時間経過に従って増加する
   請求項１または２に記載の計算装置。
8. 前記単位時間毎の算出処理において、前記更新部は、
   前記複数の要素のそれぞれについて、
   前記第１変数が前記第１値より小さい場合、前記第１変数を、前記第１値以上、前記しきい値に増加係数を加算した値以下における、乱数により定まる値に変更し、
   前記第１変数が前記第２値より大きい場合、前記第１変数を、前記しきい値に前記増加係数を加算した値以上、前記第２値以下における、乱数により定まる値に変更し、
   前記増加係数は、前記初期時刻において０であり、前記初期時刻から前記終了時刻まで時間経過に従って増加する
   請求項１または２に記載の計算装置。
9. 前記増加係数は、前記初期時刻で０、前記終了時刻で前記第２値となる時刻を変数とする一次関数、または、前記一次関数の平方根である
   請求項７または８に記載の計算装置。
10. 前記単位時間毎の算出処理において、前記更新部は、
    前記複数の要素のそれぞれについて、前記第１変数が前記第２値より大きい場合または前記第１変数が前記第１値より小さい場合、前記第２変数を、前記第２変数に乱数を乗じた値、０、予め定められた値または乱数に応じた値に変更する
    請求項１から９の何れか１項に記載の計算装置。
11. 前記単位時間毎の算出処理において、前記更新部は、
    前記複数の要素のそれぞれについて、対象時刻における前記第１変数を、前記対象時刻より単位時間前の直前時刻における前記第１変数に、前記第２変数と予め定められた定数と前記単位時間とを乗算した値を加算することにより、算出する
    請求項１から１０の何れか１項に記載の計算装置。
12. 前記単位時間毎の算出処理において、前記更新部は、
    前記複数の要素のそれぞれについて、
    前記複数の要素のそれぞれの前記第１変数と、対象要素と前記複数の要素のそれぞれとの組毎に予め定められる作用係数とに基づき、外力を算出し、
    前記直前時刻における前記第２変数に、前記外力と前記第１変数とにより定まる値に前記単位時間とを乗算した値を加算することにより、前記対象時刻における前記第２変数を算出する
    請求項１１に記載の計算装置。
13. 前記単位時間毎の算出処理において、前記更新部は、
    前記複数の要素のそれぞれについて、
    前記直前時刻における前記第１変数に、前記直前時刻における前記第２変数と前記定数と前記単位時間とを乗算した値を加算することにより、前記対象時刻における前記第１変数を算出し、
    前記複数の要素のそれぞれの前記対象時刻における前記第１変数と、前記作用係数とに基づき、前記外力を算出し、
    前記直前時刻における前記第２変数に、前記外力と前記対象時刻における前記第１変数とにより定まる値に前記単位時間とを乗算した値を加算することにより、前記対象時刻における前記第２変数を算出する
    請求項１２に記載の計算装置。
14. 前記単位時間毎の算出処理において、前記更新部は、
    前記複数の要素のそれぞれについて、
    前記複数の要素のそれぞれの前記直前時刻における前記第１変数と、前記作用係数とに基づき、前記外力を算出し、
    前記直前時刻における前記第２変数に、前記外力と前記直前時刻における前記第１変数とにより定まる値に前記単位時間とを乗算した値を加算することにより、前記対象時刻における前記第２変数を算出し、
    前記直前時刻における前記第１変数に、前記対象時刻における前記第２変数と前記定数と前記単位時間とを乗算した値を加算することにより、前記対象時刻における前記第１変数を算出する
    請求項１２に記載の計算装置。
15. 前記組合せ最適化問題は、Ｎ個の離散変数を含み、
    前記更新部は、前記Ｎ個の離散変数のうちのｉ番目の離散変数に対応するｉ番目の要素の、前記対象時刻における前記第１変数を式（１０１）または（１０２）により算出し、
    

    

    Ｎは、２以上の整数であり、
    ｉは、１からＮまでの任意の整数であり、
    Ｄは、前記定数であり、
    Δｔは、前記単位時間であり、
    ｔは、前記直前時刻であり、
    ｔ＋Δｔは、前記対象時刻であり、
    ｘ_ｉ（ｔ）は、前記ｉ番目の要素の前記直前時刻における前記第１変数であり、
    ｙ_ｉ（ｔ）は、前記ｉ番目の要素の前記直前時刻における前記第２変数であり、
    ｘ_ｉ（ｔ＋Δｔ）は、前記ｉ番目の要素の前記対象時刻における前記第１変数であり、
    ｙ_ｉ（ｔ＋Δｔ）は、前記ｉ番目の要素の前記対象時刻における前記第２変数である
    請求項１２に記載の計算装置。
16. 前記組合せ最適化問題は、ＱＵＢＯ（Ｑｕａｄｒａｔｉｃ Ｕｎｃｏｎｓｔｒａｉｎｅｄ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）問題であり、
    前記更新部は、前記ｉ番目の要素の前記対象時刻における前記第２変数を式（１０３）または式（１０４）により算出し、
    

    

    ｆ_ｉ（ｔ＋Δｔ）は、式（１０５）により表され、ｆ_ｉ（ｔ）は、式（１０６）により表され、
    

    

    ｚ_ｉ（ｔ＋Δｔ）は、式（１０７）により表され、
    ｚ_ｉ（ｔ）は、式（１０８）により表され、
    

    

    ｊは、１からＮまでの任意の整数であり、
    ｈ_ｉは、Ｎ個の局所磁場係数を含む予め定められた配列に含まれる、ｉ番目の局所磁場係数であり、
    Ｊ_ｉ，ｊは、Ｎ×Ｎ個の結合係数を含む予め定められた行列に含まれる、ｉ行、ｊ列の結合係数であり、
    ｃは、係数であり、
    ｘ_ｊ（ｔ）は、前記Ｎ個の離散変数のうちのｊ番目の離散変数に対応するｊ番目の要素の前記直前時刻における前記第１変数であり、
    ｘ_ｊ（ｔ＋Δｔ）は、前記ｊ番目の要素の前記対象時刻における前記第１変数であり、
    ｐ（ｔ）は、ｔを変数とする予め定められた関数であり、ｔが増大するに従って増加し、ｔが前記初期時刻において０となり、ｔが前記終了時刻において１となり、
    α（ｔ）は、ｔを変数とする予め定められた関数である
    請求項１５に記載の計算装置。
17. ｚ_ｉ（ｔ＋Δｔ）は、式（１０７）に代えて、式（１０９）により表され、
    ｚ_ｉ（ｔ）は、式（１０８）に代えて、式（１１０）により表され、
    

    

    ｓ_ｊ（ｔ＋Δｔ）は、前記ｊ番目の要素の前記対象時刻における前記第１変数の符号であり、
    ｓ_ｊ（ｔ）は、前記ｊ番目の要素の前記直前時刻における前記第１変数の符号である
    請求項１６に記載の計算装置。
18. ｆ_ｉ（ｔ＋Δｔ）は、式（１０５）に代えて、式（１１１）により表され、
    ｆ_ｉ（ｔ）は、式（１０６）に代えて、式（１１２）により表され、
    

    

    ｇ（ｔ＋Δｔ）は、式（１１３）により表され、ｇ（ｔ）は、式（１１４）により表され、
    

    

    ｓｇｎ（ｚ_ｉ（ｔ＋Δｔ））は、ｚ_ｉ（ｔ＋Δｔ）の符号を表し、
    ｓｇｎ（ｚ_ｉ（ｔ））は、ｚ_ｉ（ｔ）の符号を表し、
    ｚ_ｉ（ｔ＋Δｔ）は、式（１０７）に代えて、式（１１５）により表され、
    ｚ_ｉ（ｔ）は、式（１０８）に代えて、式（１１６）により表される
    

    

    請求項１６に記載の計算装置。
19. 前記組合せ最適化問題は、ＨＯＢＯ（Ｈｉｇｈｅｒ Ｏｒｄｅｒ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）問題であり、
    ｚ_ｉ（ｔ＋Δｔ）は、式（１０７）に代えて、式（１１７）により表され、
    ｚ_ｉ（ｔ）は、式（１０８）に代えて、式（１１８）により表され、
    

    

    ｋは、１からＮまでの任意の整数であり、
    Ｊ^（１） _ｉは、１階のテンソルにおける１番目の成分のｉ番目の要素であり、
    Ｊ^（２） _ｉ，ｊは、２階のテンソルにおける１番目の成分のｉ番目および２番目の成分のｊ番目の、要素であり、
    Ｊ^（３） _{ｉ，ｊ，ｋ}は、３階のテンソルにおける１番目の成分のｉ番目、２番目の成分のｊ番目および３番目の成分のｋ番目の、要素であり、
    ｘ_ｋ（ｔ）は、前記Ｎ個の離散変数のうちのｋ番目の離散変数に対応するｋ番目の要素の前記直前時刻における前記第１変数であり、
    ｘ_ｋ（ｔ＋Δｔ）は、前記ｋ番目の要素の前記対象時刻における前記第１変数である
    請求項１６または１８に記載の計算装置。
20. 前記組合せ最適化問題は、ＨＯＢＯ（Ｈｉｇｈｅｒ Ｏｒｄｅｒ Ｂｉｎａｒｙ Ｏｐｔｉｍｉｚａｔｉｏｎ）問題であり、
    ｚ_ｉ（ｔ＋Δｔ）は、式（１０９）に代えて、式（１１９）により表され、
    ｚ_ｉ（ｔ）は、式（１１０）に代えて、式（１２０）により表され、
    

    

    ｋは、１からＮまでの任意の整数を表し、
    Ｊ^（１） _ｉは、１階のテンソルにおける１番目の成分のｉ番目の要素であり、
    Ｊ^（２） _ｉ，ｊは、２階のテンソルにおける１番目の成分のｉ番目および２番目の成分のｊ番目の、要素であり、
    Ｊ^（３） _{ｉ，ｊ，ｋ}は、３階のテンソルにおける１番目の成分のｉ番目、２番目の成分のｊ番目および３番目の成分のｋ番目の、要素であり、
    ｓ_ｋ（ｔ＋Δｔ）は、前記ｋ番目の要素の前記対象時刻における前記第１変数の符号であり、
    ｓ_ｋ（ｔ）は、前記ｋ番目の要素の前記直前時刻における前記第１変数の符号である
    請求項１７に記載の計算装置。
21. 情報処理装置により、組合せ最適化問題を解くための計算方法であって、
    前記情報処理装置により、第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎の前記第１変数および前記第２変数を、前記単位時間毎に順次に且つ前記第１変数および前記第２変数を交互に更新する更新ステップと、
    前記情報処理装置により、前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力ステップと、
    を実行し、
    前記複数の要素は、前記組合せ最適化問題を表す複数の離散変数に対応し、
    前記第１変数および前記第２変数のそれぞれは、実数により表され、
    前記更新ステップにおける前記単位時間毎の算出処理において、前記情報処理装置は、前記複数の要素のそれぞれについて、
    前記第１変数を前記第２変数に基づき更新し、
    前記第２変数を前記複数の要素のそれぞれの前記第１変数に基づき更新し、
    前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値以上、予め定められたしきい値以下の値に変更し、
    前記第１変数が予め定められた第２値より大きい場合、前記第１変数を前記しきい値以上、前記第２値以下の値に変更し、
    前記第２値は、前記第１値より大きく、
    前記しきい値は、前記第１値より大きく前記第２値より小さい
    計算方法。
22. 情報処理装置を、組合せ最適化問題を解く計算装置として機能させるためのプログラムであって、
    前記情報処理装置を、
    第１変数および第２変数が対応付けられた複数の要素のそれぞれについて、初期時刻から終了時刻まで単位時間毎の前記第１変数および前記第２変数を、前記単位時間毎に順次に且つ前記第１変数および前記第２変数を交互に更新する更新部と、
    前記終了時刻における前記複数の要素のそれぞれの前記第１変数に基づき前記組合せ最適化問題の解を出力する出力部と、
    して機能させ、
    前記複数の要素は、前記組合せ最適化問題を表す複数の離散変数に対応し、
    前記第１変数および前記第２変数のそれぞれは、実数により表され、
    前記単位時間毎の算出処理において、前記更新部は、前記複数の要素のそれぞれについて、
    前記第１変数を前記第２変数に基づき更新し、
    前記第２変数を前記複数の要素のそれぞれの前記第１変数に基づき更新し、
    前記第１変数が予め定められた第１値より小さい場合、前記第１変数を、前記第１値以上、予め定められたしきい値以下の値に変更し、
    前記第１変数が予め定められた第２値より大きい場合、前記第１変数を前記しきい値以上、前記第２値以下の値に変更し、
    前記第２値は、前記第１値より大きく、
    前記しきい値は、前記第１値より大きく前記第２値より小さい
    プログラム。

Images