(57)【要約】【課題】大規模なイジング問題を高速に解くことができる計算装置、表示装置およびプログラムを提供する。【解決手段】実施形態の計算装置は、第１変数更新と第２変数更新とを含む処理手順を繰り返す処理部と、表示装置に画像を表示させる表示制御部と、を備える。処理部は、処理手順を繰り返している間は更新された第１変数ｘｉおよび第２変数ｙｉの値を出力し、処理手順の繰り返しが終了すると、最後に更新された第１変数ｘｉの符号を出力する。表示制御部は、処理部が出力する第１変数ｘｉおよび第２変数ｙｉの値に基づいて、第１変数ｘｉおよび第２変数ｙｉの少なくとも一部の時間変化を表すモニタリング画像を表示装置に表示させる。【選択図】図２
2020154524.tif 000002
【特許請求の範囲】
【請求項１】
第１変数更新と第２変数更新とを含む処理手順を繰り返す処理部と、
表示装置に画像を表示させる表示制御部と、を備え、
前記第１変数更新は、Ｎ個（Ｎは２以上の整数）の第１変数ｘｉ（ｉは１以上Ｎ以下の整数）の各々について、更新前の第１変数ｘｉに第１関数より求まる値を加えて第１変数ｘｉを更新することを含み、
前記第２変数更新は、Ｎ個の第１変数ｘｉに対応するＮ個の第２変数ｙｉの各々について、更新前の第２変数ｙｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙｉを更新することを含み、
前記第１関数の変数は、第２変数ｙｉを含み、
前記第２関数の変数は、第１変数ｘｉを含み、
前記第３関数は、行列Jｉ，ｊと第１変数ｘｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含み、
前記処理部は、前記処理手順を繰り返している間は更新された第１変数ｘｉおよび第２変数ｙｉの値を出力し、前記処理手順の繰り返しが終了すると、最後に更新された第１変数ｘｉの符号を出力し、
前記表示制御部は、前記処理部が出力する第１変数ｘｉおよび第２変数ｙｉの値に基づいて、第１変数ｘｉおよび第２変数ｙｉの少なくとも一部の時間変化を表すモニタリング画像を前記表示装置に表示させること
を特徴とする計算装置。
【請求項２】
前記表示制御部は、前記処理部が出力する第１変数ｘｉおよび第２変数ｙｉの値により表される座標点（ｘｉ，ｙｉ）の少なくとも一部を２次元座標上にプロットした前記モニタリング画像を前記表示装置に表示させること
を特徴とする請求項１に記載の計算装置。
【請求項３】
前記表示制御部は、前記処理部が出力する第１変数ｘｉおよび第２変数ｙｉの値により表される座標点（ｘｉ，ｙｉ）を、複数の小ブロックにブロック分割された２次元座標上にプロットし、各小ブロックをブロック内プロット数に応じた濃度で表した前記モニタリング画像を前記表示装置に表示させること
を特徴とする請求項１に記載の計算装置。
【請求項４】
前記表示制御部は、前記プロットされた前記座標点の集合を軌跡として表すこと
を特徴とする請求項２または３に記載の計算装置。
【請求項５】
前記表示制御部は、前記プロットされた前記座標点の間を線分で繋いで軌跡として表すこと
を特徴とする請求項４に記載の計算装置。
【請求項６】
前記表示制御部は、前記処理部が出力する第１変数ｘｉの少なくとも一部の値の時間変化を表す第１グラフと前記処理部が出力する第２変数ｙｉの少なくとも一部の値の時間変化を表す第２グラフと、の少なくともいずれかを含む前記モニタリング画像を前記表示装置に表示させること
を特徴とする請求項１に記載の計算装置。
【請求項７】
計算装置に接続される表示装置であって、
前記計算装置は、
第１変数更新と第２変数更新とを含む処理手順を繰り返し、
前記第１変数更新は、Ｎ個（Ｎは２以上の整数）の第１変数ｘｉ（ｉは１以上Ｎ以下の整数）の各々について、更新前の第１変数ｘｉに第１関数より求まる値を加えて第１変数ｘｉを更新することを含み、
前記第２変数更新は、Ｎ個の第１変数ｘｉに対応するＮ個の第２変数ｙｉの各々について、更新前の第２変数ｙｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙｉを更新することを含み、
前記第１関数の変数は、第２変数ｙｉを含み、
前記第２関数の変数は、第１変数ｘｉを含み、
前記第３関数は、行列Jｉ，ｊと第１変数ｘｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含み、
前記処理手順を繰り返している間は更新された第１変数ｘｉおよび第２変数ｙｉの値を出力し、前記処理手順の繰り返しが終了すると、最後に更新された第１変数ｘｉの符号を出力し、
前記計算装置が出力する第１変数ｘｉおよび第２変数ｙｉの値に基づいて、第１変数ｘｉおよび第２変数ｙｉの少なくとも一部の時間変化を表すモニタリング画像を表示すること
を特徴とする表示装置。
【請求項８】
コンピュータに、
第１変数更新と第２変数更新とを含む処理手順を繰り返す処理部の機能と、
表示装置に画像を表示させる表示制御部の機能と、を実現させるためのプログラムであって、
前記第１変数更新は、Ｎ個（Ｎは２以上の整数）の第１変数ｘｉ（ｉは１以上Ｎ以下の整数）の各々について、更新前の第１変数ｘｉに第１関数より求まる値を加えて第１変数ｘｉを更新することを含み、
前記第２変数更新は、Ｎ個の第１変数ｘｉに対応するＮ個の第２変数ｙｉの各々について、更新前の第２変数ｙｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙｉを更新することを含み、
前記第１関数の変数は、第２変数ｙｉを含み、
前記第２関数の変数は、第１変数ｘｉを含み、
前記第３関数は、行列Jｉ，ｊと第１変数ｘｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含み、
前記処理部は、前記処理手順を繰り返している間は更新された第１変数ｘｉおよび第２変数ｙｉの値を出力し、前記処理手順の繰り返しが終了すると、最後に更新された第１変数ｘｉの符号を出力し、
前記表示制御部は、前記処理部が出力する第１変数ｘｉおよび第２変数ｙｉの値に基づいて、第１変数ｘｉおよび第２変数ｙｉの少なくとも一部の時間変化を表すモニタリング画像を前記表示装置に表示させること
を特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、計算装置、表示装置およびプログラムに関する。
【背景技術】
【０００２】
組み合わせ最適化問題の多くは、統計力学のイジングモデルのエネルギー最小状態（基底状態）を求める問題（以下、「イジング問題」と呼ぶ）に変換できる。近年、イジング問題を高速に解く量子コンピュータ（以下、「イジングマシン」と呼ぶ）の開発が盛んに進められている。しかし、現在提案されているイジングマシンはスピン数（量子ビットの数）に限界があり、大規模なイジング問題を解くことが難しいといった課題がある。
【０００３】
一方、計算機によるシミュレーションでイジング問題を解く方法として、シミュレーテッド・アニーリングが知られている。シミュレーテッド・アニーリングは、大規模なイジング問題を扱うことが可能であるが、逐次更新アルゴリズムを採用するため並列計算に不向きであり、大規模なイジング問題を高速に解くことが難しい。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１７−７３１０６号公報
【非特許文献】
【０００５】
【非特許文献１】T．Inagaki  et  al.，Science  354，603  (2016).
【非特許文献２】H．Goto，Sci．Rep．6，21686  (2016).
【非特許文献３】Y．Haribara  et  al.，Quantum  Sci．Technol．2，044002  (2017).
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明が解決しようとする課題は、大規模なイジング問題を高速に解くことができる計算装置、表示装置およびプログラムを提供することである。
【課題を解決するための手段】
【０００７】
実施形態の計算装置は、処理部と、表示制御部と、を備える。処理部は、第１変数更新と第２変数更新とを含む処理手順を繰り返す。表示制御部は、表示装置に画像を表示させる。前記第１変数更新は、Ｎ個（Ｎは２以上の整数）の第１変数ｘｉ（ｉは１以上Ｎ以下の整数）の各々について、更新前の第１変数ｘｉに第１関数より求まる値を加えて第１変数ｘｉを更新することを含む。前記第２変数更新は、Ｎ個の第１変数ｘｉに対応するＮ個の第２変数ｙｉの各々について、更新前の第２変数ｙｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙｉを更新することを含む。前記第１関数の変数は、第２変数ｙｉを含む。前記第２関数の変数は、第１変数ｘｉを含む。前記第３関数は、行列Jｉ，ｊと第１変数ｘｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含む。前記処理部は、前記処理手順を繰り返している間は更新された第１変数ｘｉおよび第２変数ｙｉの値を出力し、前記処理手順の繰り返しが終了すると、最後に更新された第１変数ｘｉの符号を出力する。前記表示制御部は、前記処理部が出力する第１変数ｘｉおよび第２変数ｙｉの値に基づいて、第１変数ｘｉおよび第２変数ｙｉの少なくとも一部の時間変化を表すモニタリング画像を前記表示装置に表示させる。
【図面の簡単な説明】
【０００８】
【図１】図１は、実施形態に係る計算装置のハードウェア構成例を示すブロック図である。
【図２】図２は、実施形態に係る計算装置の機能的な構成例を示す図である。
【図３】図３は、実施形態に係る計算装置の動作の一例を示すフローチャートである。
【図４】図４は、処理部による処理の一例を示すフローチャートである。
【図５】図５は、処理部による処理の一例を示すフローチャートである。
【図６】図６は、処理部による処理の一例を示すフローチャートである。
【図７】図７は、処理部による処理の一例を示すフローチャートである。
【図８】図８は、処理部による処理の一例を示すフローチャートである。
【図９】図９は、処理部による処理の一例を示すフローチャートである。
【図１０】図１０は、処理部による処理の一例を示すフローチャートである。
【図１１】図１１は、処理部による処理の一例を示すフローチャートである。
【図１２】図１２は、モニタリング画像の一例を示す図である。
【図１３】図１３は、モニタリング画像の一例を示す図である。
【図１４】図１４は、モニタリング画像の一例を示す図である。
【図１５】図１５は、モニタリング画像の一例を示す図である。
【図１６】図１６は、モニタリング画像の一例を示す図である。
【図１７】図１７は、モニタリング画像の一例を示す図である。
【発明を実施するための形態】
【０００９】
以下、本発明を適用した計算装置、表示装置およびプログラムの具体的な実施形態について、添付図面を参照しながら詳細に説明する。なお、以下の説明では、共通の要素については同一の参照符号を付して、重複した説明を適宜省略する。
【００１０】
図１は、本実施形態に係る計算装置１００のハードウェア構成例を示すブロック図である。計算装置１００は、図１に示すように、第１プロセッサ回路１０１と、第２プロセッサ回路１０２と、主記憶装置１０３と、補助記憶装置１０４と、機器Ｉ／Ｆ１０５と、通信Ｉ／Ｆ１０６とを備える。
【００１１】
第１プロセッサ回路１０１は、例えばＧＰＵ（Graphics  Processing  Unit）を含む。第１プロセッサ回路１０１は、計算量が多い特定の処理を並列計算によって高速に行う。本実施形態において、特定の処理とは、後述の新規な運動方程式の時間発展を計算することによってイジング問題を解く（近似解を求める）処理である。
【００１２】
第２プロセッサ回路１０２は、例えばＣＰＵ（Central  Processing  Unit）を含む。第２プロセッサ回路１０２は、計算装置１００を適切に動作させるための制御など、比較的複雑な処理を行う。本実施形態では、例えば、ユーザが入力した各種パラメータの設定や変数の設定、後述のモニタリング画像の表示制御などが、第２プロセッサ回路１０２により行われる。
【００１３】
主記憶装置１０３は、例えばＲＡＭ（Random  Access  Memory）を含む。主記憶装置１０３は、第１プロセッサ回路１０１や第２プロセッサ回路１０２が作業領域として利用する記憶装置である。
【００１４】
補助記憶装置１０４は、例えばＨＤＤ（Hard  Disk  Drive）あるいはＳＳＤ（Solid  State  Drive）などの大容量ストレージ装置を含む。補助記憶装置１０４は、計算装置１００において使用されるプログラムやデータなどを格納する。補助記憶装置１０４に格納されたプログラムやデータは、第１プロセッサ回路１０１や第２プロセッサ回路１０２によって読み出されて主記憶装置１０３上にロードされる。第１プロセッサ回路１０１や第２プロセッサ回路１０２は、例えば主記憶装置１０３上にロードしたプログラムに従って、上述の処理を実行することができる。
【００１５】
機器Ｉ／Ｆ１０５は、計算装置１００に対して表示装置１１０（例えば液晶ディスプレイ）や入力装置１２０（例えばキーボード、マウス、タッチ式入力パネル、音声認識入力装置など）などの外部機器を接続するためのインタフェースである。
【００１６】
通信Ｉ／Ｆ１０６は、計算装置１００をネットワーク１３０に接続するためのインタフェースである。ネットワーク１３０は、例えばインターネットやＬＡＮ（Local  Area  Network）などを含む。
【００１７】
本実施形態に係る計算装置１００としては、例えば汎用のコンピュータを用いることができる。また、互いに接続された複数のコンピュータを用いて計算装置１００を構成してもよい。また、計算装置１００の少なくとも一部（例えば第１プロセッサ回路１０１）として、例えばＡＳＩＣ（Application  Specific  Integrated  Circuit）やＦＰＧＡ（Field-Programmable  Gate  Array）などの専用の回路を用いてもよい。
【００１８】
図２は、本実施形態に係る計算装置１００の機能的な構成例を示すブロック図である。計算装置１００は、図２に示すように、機能的な構成要素として、設定部１０と、処理部２０と、表示制御部３０とを備える。設定部１０および表示制御部３０は、例えば第２プロセッサ回路１０２により実現することができる。処理部２０は、例えば第１プロセッサ回路１０１により実現することができる。
【００１９】
設定部１０は、後述の新規な運動方程式を解くための各種パラメータの設定や変数の初期化などを行う。各種パラメータは、例えば、ユーザが入力装置１２０を用いて入力した値に設定される。変数は、イジングモデルのスピン数に対応するＮ個の（Ｎは２以上の整数）の第１変数ｘｉ（ｉは１以上Ｎ以下の整数）と、Ｎ個の第１変数ｘｉに対応するＮ個の第２変数ｙｉとを含む。第１変数ｘｉは、イジングスピンの値域を実数に拡張した変数である。第２変数ｙｉは、後述の新規な運動方程式中に現れる値域が実数の変数である。Ｎ個の第１変数ｘｉの集合を第１変数群｛ｘ｝、Ｎ個の第２変数ｙｉの集合を第２変数群｛ｙ｝と呼ぶ。これら第１変数ｘｉおよび第２変数ｙｉは、計算を開始する前に、例えば絶対値が０．１以下の乱数によってランダムに初期化される。
【００２０】
処理部２０は、後述の新規な運動方程式の時間発展を計算するために、第１変数更新と第２変数更新とを含む処理手順を繰り返す。第１変数更新は、Ｎ個の第１変数ｘｉの各々について、更新前の第１変数ｘｉに第１関数より求まる値を加えて第１変数ｘｉを更新することを含む。第２変数更新は、Ｎ個の第２変数ｙｉの各々について、更新前の第２変数ｙｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙｉを更新することを含む。第１関数の変数は、第２変数ｙｉを含み、第２関数の変数は、第１変数ｘｉを含む。第３関数は、行列Jｉ，ｊと第１変数ｘｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含む。第１関数、第２関数および第３関数は、後述の新規な運動方程式から導かれる関数である。これらの具体例については詳細を後述する。
【００２１】
処理部２０は、上述の処理手順を繰り返している間は、更新された第１変数ｘｉおよび第２変数ｙｉの値を出力する。そして、所定の終了条件を満たし、上述の処理手順の繰り返しが終了すると、処理部２０は、最後に更新された第１変数ｘｉの符号「±１」を出力する。後述の新規な運動方程式の時間発展が適切に計算された場合、最後に更新されたＮ個の第１変数ｘｉの符号「±１」は、基底状態のイジングスピンを表しているものと期待される。
【００２２】
表示制御部３０は、処理部２０が出力する第１変数ｘｉおよび第２変数ｙｉの値に基づいて、第１変数ｘｉおよび第２変数ｙｉの少なくとも一部の時間変化を表すモニタリング画像を表示装置１１０に表示させる。モニタリング画像の具体例については詳細を後述する。
【００２３】
本実施形態に係る計算装置１００によるシミュレーションでイジング問題を解く場合、処理部２０による処理手順の繰り返しによって後述の新規な運動方程式の時間発展が進むにつれて、Ｎ個の第１変数ｘｉの値が０付近からプラス・マイナスの一定値に分岐して収束していくことが理想である。しかし、パラメータの設定が適切でない場合は、Ｎ個の第１変数ｘｉの値がうまく分岐しなかったり無限大に発散したりする場合がある。
【００２４】
そこで、本実施形態では、表示制御部３０が、処理手順を繰り返している間に処理部２０から出力される第１変数ｘｉおよび第２変数ｙｉの値に基づいて、第１変数ｘｉおよび第２変数ｙｉの少なくとも一部の時間変化を表すモニタリング画像（新規な運動方程式の時間発展に伴う第１変数ｘｉおよび第２変数ｙｉの時間変化を可視化した画像）を生成する。そして、表示制御部３０は、このモニタリング画像を表示装置１１０に表示させるようにしている。したがって、ユーザは、表示装置１１０に表示されているモニタリング画像を参照することにより、第１変数ｘｉの分岐が適切に行われているか否かを確認することができ、パラメータの設定が適切であるかが確認できる。また、分岐が適切でない場合に、モニタリング画像で表される第１変数ｘｉおよび第２変数ｙｉの時間変化の様子から、パラメータのチューニングの方向性を判断することができる。
【００２５】
次に、本実施形態に係る計算装置１００の動作について説明する。図３は、計算装置１００の動作の一例を示すフローチャートである。
【００２６】
計算装置１００の動作が開始されると、まず、設定部１０が、後述の新規な運動方程式を解くための各種パラメータの設定を行う（ステップＳ１）。パラメータは、例えば、第１変数ｘｉや第２変数ｙｉの数を表す「Ｎ」、後述の「Ｔ」、「Ｄ」、「Ｋ」、「ｃ」、「ｐ」などを含む。また、パラメータは、要素間の相互作用の強さを表す行列「Jｉ，ｊ」（第１パラメータ）を含む。また、パラメータは、外場を表すベクトル「ｈｉ」（第２パラメータ）をさらに含んでもよい。
【００２７】
次に、設定部１０は、第１変数群｛ｘ｝に含まれる各第１変数ｘｉおよび第２変数群｛ｙ｝に含まれる各第２変数ｙｉを、例えば絶対値が０．１以下の乱数によってランダムに初期化する（ステップＳ２）。
【００２８】
次に、処理部２０が、後述の新規な運動方程式の時間発展の計算を行う（ステップＳ３）。すなわち、処理部２０は、上述の第１変数更新と第２変数更新とを含む処理手順を複数回繰り返すことにより、後述の新規な運動方程式の時間発展を計算する。第１変数更新では、第２変数群｛ｙ｝を用いて第１変数群｛ｘ｝が更新され、第２変数更新では、第１変数群｛ｘ｝を用いて第２変数群｛ｙ｝が更新される。複数回の処理手順の少なくとも一部において、第１変数更新の後に第２変数更新が行われてもよいし、第２変数更新の後に第１変数更新が行われてもよい。第１変数更新と第２変数更新とを含む処理手順は、所定の終了条件が満たされるまで繰り返される。
【００２９】
処理手順の繰り返しが終了すると、処理部２０は、最後に更新された第１変数群｛ｘ｝に含まれる各第１変数ｘｉの符号を出力する（ステップＳ４）。これら第１変数ｘｉの符号は、基底状態のイジングスピンを表している。本実施形態に係る計算装置１００は、以上のように、新規な運動方程式の時間発展を計算することにより、大規模なイジング問題を高速に解くことができる。
【００３０】
ここで、本実施形態に係る計算装置１００が扱うイジング問題について説明する。イジングエネルギーは、例えば下記式（１）のように表される。
【数１】
2020154524.tif 000003
この式（１）において、「Ｎ」はイジングスピンの数（スピン数）である。「ｓｉ」はｉ番目のイジングスピンであり、例えば「ｓｉ」＝±１である。「Jｉ，ｊ」はイジングスピン間の相互作用の強さを表す行列（第１パラメータ）であり、例えば、対角成分（対角要素）が全てゼロの実対象行列である。「ｈｉ」は外場を表すベクトル（第２パラメータ）である。
【００３１】
上記式（１）に関して、量子分岐マシンの古典モデル（以下、「古典分岐マシン」と呼ぶ）が提案されている。例えば特許文献１には、古典分岐マシンによるシミュレーションの一例が記載されている。古典分岐マシンにおいては、運動方程式は、以下の式（２）乃至式（４）の連立常微分方程式で与えられる。
【数２】
2020154524.tif 000004
【数３】
2020154524.tif 000005
【数４】
2020154524.tif 000006
【００３２】
これら式（２）乃至式（４）において、「Ｎ」は、例えば、イジングスピンの数（スピン数）に対応する。「Ｄ」は、例えば、特許文献１に記載の「離調」に対応するパラメータである。「ｃ」は定数である。「ｐ」は、例えば、特許文献１に記載の「ポンプ振幅（ポンプレート）」に対応するパラメータである。「Ｋ」は、例えば、特許文献１に記載の「カー係数」に対応するパラメータである。上記式（２）乃至式（４）で表される運動方程式において、第２パラメータ「ｈｉ」は設けられなくてもよい。その場合は、上記式（３）および式（４）の中の「ｈｉ」を含む項は無視される。
【００３３】
上記式（２）乃至式（４）において、ｐ（ｔ）をゼロから十分大きな値へ増加させたときの「ｘｉ」の最終値（最後に更新された値）の符号「±１」が、最適解（基底状態）のイジングスピン「ｓｉ」となる。「ａ（ｔ）」は、「ｐ（ｔ）」とともに増加するパラメータである。「ａ（ｔ）」は、例えば、下記式（５）により表される。
【数５】
2020154524.tif 000007
【００３４】
上記の古典分岐マシンは、上記式（２）乃至式（４）において、「Ｈ」をハミルトニアンとしたハミルトン力学系であると見なすことができる。
【００３５】
計算機によるシミュレーションでイジング問題を解く方法としては、シミュレーテッド・アニーリングが広く知られている。この方法では、逐次更新アルゴリズムが採用される。逐次更新アルゴリズムにおいては、複数のスピンが１つずつ更新される。このような逐次更新アルゴリズムは、並列計算に向かない。
【００３６】
これに対して、上記の古典分岐マシンの運動方程式をデジタル計算機で離散解法によって解くことが考えられる。このアルゴリズムは、シミュレーテッド・アニーリングとは異なり、並列更新アルゴリズムである。並列更新アルゴリズムにおいては、複数の変数を同時に更新することができる。このため、並列計算による高速化が期待できる。
【００３７】
上記式（２）乃至式（４）を用いるアプローチには、以下の課題があると考えられる。すなわち、最も計算量が大きい行列Jｉ，ｊに関する積和計算が、第１変数ｘｉの更新と第２変数ｙｉの更新の双方において必要となる。上記の運動方程式は、数値的に容易には解けないため、例えば、計算量が大きい離散解法（例えば４次のルンゲ・クッタ法など）が必要となる。このため、並列計算を導入しても十分な高速化が図れない懸念がある。
【００３８】
そこで、本実施形態においては、上記式（２）乃至式（４）で表される運動方程式ではなく、例えば、下記式（６）乃至式（８）で表される新規な運動方程式を解く。
【数６】
2020154524.tif 000008
【数７】
2020154524.tif 000009
【数８】
2020154524.tif 000010
【００３９】
これら式（６）乃至式（８）において、「Ｎ」は、例えば、イジングスピンの数（スピン数）に対応する。「Ｄ」は、例えば、特許文献１に記載の「離調」に対応するパラメータである。「ｃ」は定数である。「ｐ」は、例えば、特許文献１に記載の「ポンプ振幅（ポンプレート）」に対応するパラメータ（演算パラメータ）である。「Ｋ」は、例えば、特許文献１に記載の「カー係数」に対応するパラメータである。これらのパラメータは、例えば、図３のステップＳ１において設定される。１つの例において、「Ｋ」は１に設定される。上記式（６）乃至式（８）で表される運動方程式において、第２パラメータ「ｈｉ」は設けられなくてもよい。その場合は、上記式（７）および式（８）の中の「ｈｉ」を含む項は無視される。
【００４０】
上記式（６）乃至式（８）で表される新規な運動方程式を用いる場合、最も計算量が大きい行列Jｉ，ｊに関する積和演算は、第２変数ｙｉの更新においてのみ行われ、第１変数ｘｉの更新では行われない。したがって、計算量が削減される。この新規な運動方程式では、第１変数ｘｉの時間微分から「ｐ」が消去されている。この新規な運動方程式において、第１変数ｘｉの時間微分は、第２変数ｙｉを含み、第１変数ｘｉを含まない。また、第２変数ｙｉの時間微分は、第１変数ｘｉを含み、第２変数ｙｉを含まない。ハミルトニアンにおいて、第１変数ｘｉおよび第２変数ｙｉは、互いに分離されている。このため、計算量が小さく安定な離散解法が適用可能となる。例えば、シンプレクティック・オイラー法と呼ばれる方法が適用できる。
【００４１】
イジング問題を解くアプローチとして新規な運動方程式を用いた場合に、高い性能（例えば、高い精度）が維持できることが分かった。本実施形態に係る計算装置１００では、上記の分離可能なハミルトニアンを持つハミルトン力学系（新しい古典分岐マシン）の新規な運動方程式が、例えば、シンプレクティック・オイラー法で解かれる。本実施形態に係る計算装置１００は、このような新しいアルゴリズムの計算が並列計算によって、できる限り高速に実行されるように構成される。
【００４２】
実施形態においては、図３のステップＳ３において、処理部２０が第１変数更新および第２変数更新を含む処理手順を繰り返すことで上述の新規な運動方程式の時間発展を計算し、イジングモデルの基底状態を求める。例えば、「ｐ（ｔ）」をゼロから十分大きな値へ増加させたときの第１変数ｘｉの最終値の符号「±１」が、基底状態のイジングスピンｓｉとなる。以下では、図３のステップＳ３における処理の詳細について説明する。
【００４３】
図４は、処理部２０による処理の一例を示すフローチャートである。この図４のフローチャートで示す処理は、図３のステップＳ３においてサブルーチンとして実行される。
【００４４】
このサブルーチンが開始されると、処理部２０は、まず、時刻「ｔ」と演算パラメータ「ｐ」および「ａ」を初期化する（ステップＳ１０１）。１つの例において、時刻「ｔ」と演算パラメータ「ｐ」および「ａ」は、０に初期化される。なお、「Ｔ」は、時刻「ｔ」の最終値に対応する。「ｄｔ」は、時刻「ｔ」の１ステップ当たりの増加量である。「ｄｐ」は、演算パラメータ「ｐ」の１ステップ当たりの増加量である。これら「Ｔ」、「ｄｔ」、「ｄｐ」は、例えば、図３のステップＳ１において設定されるパラメータに含まれる。
【００４５】
処理部２０は、「ｔ」が「Ｔ」よりも小さいときに、以下に説明する第１変数更新および第２変数更新を含む処理手順を実施する（ステップＳ１０２）。すなわち、第１変数更新および第２変数更新を含む処理手順は、「ｔ」が「Ｔ」以上という終了条件を満たすまで繰り返されるループ処理の１ステップとして実施される。なお、「ｐ」が予め適切に設定された「Ｐ」以上であることを終了条件とし、「ｐ」が「Ｐ」よりも小さいときに、第１変数更新および第２変数更新を含む処理手順を実施するようにしてもよい。
【００４６】
ループ処理の１ステップとして実施される処理手順において、処理部２０は、例えば、まず、第１変数更新を行う（ステップＳ１１０）。第１変数更新では、Ｎ個のイジングスピンｓｉに対応するＮ個の第１変数ｘｉの更新が行われる。例えば、更新前の第１変数ｘｉに、ｄｔ＊Ｄ＊ｙｉを加えて得られる値を、更新後の第１変数ｘｉとする。ここで、「＊」は、積の記号である。
【００４７】
次に、処理部２０は、第２変数更新を行う（ステップＳ１２０）。第２変数更新では、Ｎ個の第１変数ｘｉに対応するＮ個の第２変数ｙｉの更新が行われる。第２変数更新は、第１サブ更新（ステップＳ１２１）と第２サブ更新（ステップＳ１２２）とを含む。第１サブ更新と第２サブ更新の順番は、入れ替えが可能である。第１サブ更新の少なくとも一部と第２サブ更新の少なくとも一部とが同時に行われてもよい。
【００４８】
第１サブ更新では、例えば、更新前の第２変数ｙｉに、ｄｔ＊［（ｐ−Ｄ−ｘｉ＊ｘｉ）＊ｘｉ−ｃ＊ｈｉ＊ａ］を加えて得られる値を、更新後の第２変数ｙｉとする。
【００４９】
第２サブ更新では、例えば、更新前の第２変数ｙｉに、ｄｔ＊ｃ＊Σ（Ｊｉ，ｊ＊ｘｊ）を加えて得られる値を更新後の第２変数ｙｉとする。「Σ」は、行列Ｊに関する和を表す。例えば、「ｄｔ＊ｃ＊Ｊ」を行列Ｊとしてもよい。この場合、「ｄｔ＊ｃ＊」の演算を実際に行わなくてもよい。
【００５０】
第１変数更新および第２変数更新が終了すると、処理部２０は、更新後の第１変数ｘｉおよび第２変数ｙｉの値を出力する（ステップＳ１３０）。これら更新後の第１変数ｘｉおよび第２変数ｙｉの値は、表示制御部３０が後述のモニタリング画像を生成して表示装置１１０に表示させるために利用される。
【００５１】
次に、処理部２０は、時刻「ｔ」と演算パラメータ「ｐ」および「ａ」を更新する（ステップＳ１４０）。例えば、更新前の「ｔ」に「ｄｔ」を加えて得られる値を、更新後の「ｔ」とする。また、更新前の「ｐ」に「ｄｐ」を加えて得られる値を、更新後の「ｐ」とする。また、「ｐ１／２」を、更新後の「ａ」とする。本実施形態では、「ｐ」の更新において、任意の増加関数を用いた線形増加が適用される。
【００５２】
そして、更新後の「ｔ」が「Ｔ」よりも小さければ、ステップＳ１０２に戻って以降の処理手順を繰り返す（ステップＳ１０３）。その後、更新後の「ｔ」が「Ｔ」以上になると、ループ処理を終了して図３のステップＳ４に進む。なお、更新後の「ｐ」が予め適切に設定された「Ｐ」よりも小さければステップＳ１０２に戻り、更新後の「ｐ」が「Ｐ」以上になると、ループ処理を終了して図３のステップＳ４に進むようにしてもよい。
【００５３】
なお、以上の説明では、ループ処理の各ステップにおいて、第１変数ｘｉおよび第２変数ｙｉの更新が、第１変数更新（ステップＳ１１０）、第２変数更新の第１サブ更新（ステップＳ１２１）、第２変数更新の第２サブ更新（ステップＳ１２２）の順番で行われるものとしたが、これらの処理の順番は任意に変更することができる。例えば、ループ処理の少なくとも一部のステップにおいて、第１変数ｘｉおよび第２変数ｙｉの更新を、図５に示すように、第２変数更新の第１サブ更新（ステップＳ１２１）、第２変数更新の第２サブ更新（ステップＳ１２２）、第１変数更新（ステップＳ１１０）の順番で行うようにしてもよい。
【００５４】
また、ループ処理の少なくとも一部のステップにおいて、第１変数ｘｉおよび第２変数ｙｉの更新を、図６に示すように、第２変数更新の第１サブ更新（ステップＳ１２１）、第１変数更新（ステップＳ１１０）、第２変数更新の第２サブ更新（ステップＳ１２２）の順番で行うようにしてもよい。また、ループ処理の少なくとも一部のステップにおいて、第１変数ｘｉおよび第２変数ｙｉの更新を、図７に示すように、第２変数更新の第２サブ更新（ステップＳ１２２）、第１変数更新（ステップＳ１１０）、第２変数更新の第１サブ更新（ステップＳ１２１）の順番で行うようにしてもよい。
【００５５】
本実施形態に係る計算装置１００で実行されるアルゴリズムは、例えば、以下を含む。すなわち、行列Ｊｉ，ｊを取得、または、計算により定める。このとき、ベクトルｈｉをさらに取得、または、計算により定めるようにしてもよい。行列Ｊｉ，ｊおよびベクトルｈｉは、イジングモデルにおけるパラメータである。２種類の変数（第１変数群｛ｘ｝および第２変数群｛ｙ｝）を用いる。一方の変数の更新には、他方の変数の値を用いる。一方の変数の更新において、その一方の変数の値を用いない。一方の変数の更新の後に、更新後のその一方の変数の値を用いて、他方を更新する。
【００５６】
すなわち、第１変数群｛ｘ｝を構成する各第１変数ｘｉの更新（第１変数更新）は、更新前の第１変数ｘｉに第１関数より求まる値を加えて第１変数ｘｉを更新することを含む。また、第２変数群｛ｙ｝を構成する各第２変数ｙｉの更新（第２変数更新）は、更新前の第２変数ｙｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙｉを更新することを含む。第１関数は、第１変数群｛ｘ｝から独立している。第１変数群｛ｘ｝の値を変更しても、第１関数の結果の値は変化しない。第２関数は、第２変数群｛ｙ｝から独立している。第２変数群｛ｙ｝の値を変更しても、第２関数の結果の値は変化しない。第３関数は、第２変数群｛ｙ｝から独立している。第２変数群｛ｙ｝の値を変更しても、第３関数の結果の値は変化しない。
【００５７】
第１関数は、例えば、ｄｔ＊Ｄ＊ｙｉである（図４参照）。第２関数は、例えば、ｄｔ＊［（ｐ−Ｄ−ｘｉ＊ｘｉ）＊ｘｉ−ｃ＊ｈｉ＊ａ］である（図４参照）。第３関数は、例えば、ｄｔ＊ｃ＊Σ（Ｊｉ，ｊ＊ｘｊ）である（図４参照）。第１関数の変数は、第２変数ｙｉを含む。第２関数の変数は、第１変数ｘｉを含む。第３関数は、行列Ｊｉ，ｊと第１変数ｘｊとの積和演算Σ（Ｊｉ，ｊ＊ｘｊ）を含む。
【００５８】
第２関数は、例えば、第１変数ｘｉの非線形関数である第４関数を含む。第４関数は、演算パラメータ「ｐ」も含む。２種類の変数の更新とともに「ｐ」は変化する。２種類の変数の更新とともに「ｐ」が変化すると、第４関数の実数根の数が変化する。「関数の実数根」とは、関数の値がゼロとなる変数の値（実数）のことである。第２変数更新において第４関数のみを考慮した場合、第４関数の実数根は非線形力学系の固定点（fixed  point）に対応する。（ハミルトン力学系では、固定点はハミルトニアンの極値に対応する。）よって、第４関数の実数根の数が変化することは、固定点の数が変化することに対応する。これは、非線形力学系の分岐現象に対応する。
【００５９】
本実施形態に係る計算装置１００で用いるアルゴリズムでは、第１変数ｘｉおよび第２変数ｙｉの初期値が、初期の１つの安定固定点（stable  fixed  point）の近傍に設定される。「ｐ」を変化させることで分岐（bifurcation）を起こす。分岐後の複数の安定な固定点（変数の値はこの複数の安定固定点のうちの１つの近傍へと変化する）と、解きたい組合せ最適化問題の離散変数とを対応させる。これにより、分岐現象によって、その組合せ最適化問題を解く。例えば、上記の例では、分岐後の安定固定点における各第１変数ｘｉの値が、正負の２値であり、その符号と、イジングスピン（イジング問題の離散変数）とが対応付けられている。初期の安定固定点が原点であるため、各第１変数ｘｉと各第２変数ｙｉの初期値を原点の近傍（つまり、絶対値が０．１以下の小さな乱数）の値に設定する。
【００６０】
第４関数は、例えば、ｄｔ＊（ｐ−Ｄ’−ｘｉ＊ｘｉ）＊ｘｉである。「Ｄ’」は、０≦Ｄ’≦Ｄを満たす適切な定数である。初期時刻ではｐ＝０であり、第４関数の根はｘｉ＝０の１つだけであるが、ｐがＤ’よりも大きくなると、根は３つとなり、そのうちの正負２つの根が、イジングスピンと対応付けられる。なお、第２関数ｄｔ＊［（ｐ−Ｄ−ｘｉ＊ｘｉ）＊ｘｉ−ｃ＊ｈｉ＊ａ］は、ｄｔ＊（ｐ−Ｄ’−ｘｉ＊ｘｉ）＊ｘｉ＋ｄｔ＊［−（Ｄ−Ｄ’）＊ｘｉ−ｃ＊ｈｉ＊ａ］のように、第４関数と１次関数との和で表すことができる。したがって、第２関数は第４関数を含む。
【００６１】
上記の第４関数は、例えば、３次関数である。このような処理により、例えば、ニューラルネットワークで用いられている非線形関数（例えばシグモイド関数）を用いた計算よりも、計算が容易になる。
【００６２】
本実施形態に係る計算装置１００においては、時間ステップ（例えば「ｄｔ」）を大きくすると、計算が高速になる。一方、時間ステップを過度に大きくすると、計算が不安定になる。このことを考慮して、計算の一部における時間ステップを大きくし、計算の別の一部では、時間ステップを小さくてもよい。例えば、計算量の大きい、行列Ｊｉ，ｊと第１変数ｘｊとの積和演算を含む第２変数更新の第２サブ更新（ステップＳ１２２）には大きい時間ステップを適用し、その他の第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）には小さい時間ステップを適用するようにしてもよい。これにより、さらなる高速化が可能である。
【００６３】
このような計算を行う場合の処理部２０による処理の一例を図８に示す。図８に示す例では、１つのループ（ステップＳ１０２〜ステップＳ１０３）の中に、小さいループ（ステップＳ１０４〜ステップＳ１０５）が設けられている。ステップＳ１０４において、ループ変数「ｍ」は、１以上Ｍ以下である。小さいループの中において、第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）がＭ回繰り返される。これら第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）をＭ回繰り返した後、ステップＳ１２２に進んで第２変数更新の第２サブ更新が行われる。
【００６４】
なお、小さいループ中における第１変数更新（ステップＳ１１０）と第２変数更新の第１サブ更新（ステップＳ１２１）の順序は入れ替えが可能である。第１変数更新（ステップＳ１１０）と第２変数更新の第１サブ更新（ステップＳ１２１）の順序を入れ替えた例を図９に示す。
【００６５】
なお、図８および図９の例では、第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）をＭ回繰り返した後に第２変数更新の第２サブ更新（ステップＳ１２２）を行うようにしているが、図１０および図１１に示すように、第２変数更新の第２サブ更新（ステップＳ１２２）を行った後に、第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）をＭ回繰り返すようにしてもよい。
【００６６】
本実施形態に係る計算装置１００おいて、処理部２０は、主記憶装置１０３などのメモリに保持されたデータを読み出し、データを更新し、更新されたデータをメモリに保持させる。例えば、上述の第１変数更新（ステップＳ１１０）において、処理部２０は、更新前の第１変数ｘｉと第２変数ｙｉをメモリから取得し、第２変数ｙｉを用いて第１関数を計算し、得られた値を更新前の第１変数ｘｉに加えることにより第１変数ｘｉを更新し、更新された第１変数ｘｉをメモリに保持させる。
【００６７】
また、上述の第２変数更新の第１サブ更新（ステップＳ１２１）において、処理部２０は、更新前の第２変数ｙｉと第１変数ｘｉをメモリから取得し、第１変数ｘｉを用いて第２関数を計算し、得られた値を更新前の第２変数ｙｉに加えることにより第２変数ｙｉを更新し、更新された第２変数ｙｉをメモリに保持させる。また、上述の第２変数更新の第２サブ更新（ステップＳ１２２）において、処理部２０は、更新前の第２変数ｙｉと行列Ｊｉ，ｊと第１変数ｘｊをメモリから取得し、行列Ｊｉ，ｊと第１変数ｘｊを用いて第３関数を計算し、得られた値を更新前の第２変数ｙｉに加えることにより第２変数ｙｉを更新し、更新された第２変数ｙｉをメモリに保持させる。
【００６８】
実施形態において、例えば、行列Ｊｉ，ｊが疎行列の場合、疎行列の圧縮形式を用いてもよい。疎行列の圧縮形式には、例えば、ＣＯＯ（coordinate）形式またはＣＳＲ（compressed  sparse  row）形式などが適用できる。疎行列の圧縮形式を用いることで、例えば、メモリサイズを節約できる。疎行列の圧縮形式を用いることで、例えば、行列Ｊｉ，ｊと第１変数ｘｊとの積和演算を高速に実施できる。
【００６９】
以下、定数「ｃ」の例について述べる。例えば、離調「Ｄ」は、行列Ｊｉ，ｊの最大固有値λｍａｘのｃ倍よりも大きくされる（例えば、非特許文献２参照）。「Ｄ」が大きすぎると、無駄な計算時間が生じる。このため、例えば、「Ｄ」がλｍａｘのｃ倍と実質的に等しいように設定する。この場合、ｃ＝Ｄ／λｍａｘとなる。一方、１つの例において、行列Ｊｉ，ｊは、実対称行列である。この場合に、行列Ｊｉ，ｊのサイズが十分大きい場合、λｍａｘは、２σ×Ｎ１／２と実質的に同じになる。この関係は、ランダム行列のウィグナーの半円則に基づく。「σ」は、行列Ｊの非対角成分の標準偏差である。この場合、ｃ＝Ｄ／（２σ×Ｎ１／２）と設定するとよい。
【００７０】
本実施形態において、精度を向上させる方法として、上記式（７）の非線形関数として用いられる関数を変更することが考えられる。例えば、上記式（７）の代わりに下記式（９）を用いてもよい。
【数９】
2020154524.tif 000011
この式（９）において、「ｎ」は、２以上の偶数である。このような関数を用いることで、例えば、イジング問題の解の精度を向上することができる。
【００７１】
次に、本実施形態に係る計算装置１００によるモニタリング画像の表示制御について説明する。モニタリング画像は、処理部２０が上述の新規な運動方程式の時間発展を計算するために処理手順を繰り返している間の第１変数ｘｉおよび第２変数ｙｉの時間変化を可視化した画像である。モニタリング画像は、上述のループ処理の１ステップごとに処理部２０が出力する第１変数ｘｉおよび第２変数ｙｉの値に基づいて、表示制御部３０により生成され、表示装置１１０に表示される。
【００７２】
表示制御部３０は、例えば、上述のループ処理の１ステップごとに処理部２０が出力する第１変数ｘｉおよび第２変数ｙｉの値を、予め定めたステップ数ｍ分蓄積する。そして、蓄積した各ステップにおける第１変数ｘｉおよび第２変数ｙｉの値により表される座標点（ｘｉ，ｙｉ）を２次元座標上にプロットすることによりモニタリング画像を生成し、このモニタリング画像を表示装置１１０に表示させる。この際、モニタリング画像上の座標点（ｘｉ，ｙｉ）をｉの値に応じてＮ種類に色分けして表示させることが望ましい。
【００７３】
このモニタリング画像では、２次元座標上における座標点（ｘｉ，ｙｉ）のステップごとの位置変化を軌跡として表すことができ、上述の新規な運動方程式の時間発展に伴う第１変数ｘｉおよび第２変数ｙｉの時間変化を視覚的に分かり易く表示することができる。ここで、色分けしない場合であっても、それぞれのｉについて過去に蓄積されたｍ個の点を表示して、それらのｍ個の点の集合を軌跡で描くようにする、または、それらｍ個の点を順に線分で繋いで軌跡として表す、という方法を採ることができる。なお、モニタリング画像の表示中、処理部２０が新たに出力する第１変数ｘｉおよび第２変数ｙｉの値は随時蓄積され、予め定めたステップ数分蓄積されると新たにモニタリング画像が生成され、表示装置１１０が表示するモニタリング画像が更新される。モニタリング画像の更新単位となるステップ数は、上述の新規な運動方程式の時間発展の初期で第１変数ｘｉが未分岐の状態と、上述の新規な運動方程式の時間発展が進んで第１変数ｘｉが分岐していることが期待される状態とが混在しないように設定するとよい。
【００７４】
図１２乃至図１４は、モニタリング画像の一例を示す図であり、ｉ＝１〜４の４つの座標点（ｘｉ，ｙｉ）を各ステップごとに２次元座標上にプロットしたモニタリング画像を示している。図１２に示すモニタリング画像は、上述の新規な運動方程式の時間発展の初期状態に対応し、４つの座標点（ｘｉ，ｙｉ）の軌跡が２次元座標の原点（０，０）付近に集中していることから、第１変数ｘｉが未分岐の状態であることが分かる。
【００７５】
図１３に示すモニタリング画像は、上述の新規な運動方程式の時間発展が進んで、図１２に示すモニタリング画像から更新されたモニタリング画像の一例である。この図１３の例では、４つの座標点（ｘｉ，ｙｉ）の軌跡がすべて２次元座標のＸ軸のプラス側またはマイナス側の一定値付近にあり、上述の新規な運動方程式の時間発展が進んだことにより４つの第１変数ｘｉの値がうまく分岐していることが分かる。したがって、ユーザはこのモニタリング画像を参照することで、パラメータの設定が適切であることを確認できる。
【００７６】
図１４に示すモニタリング画像は、上述の新規な運動方程式の時間発展が進んで、図１２に示すモニタリング画像から更新されたモニタリング画像の他の例である。この図１４の例では、４つの座標点（ｘｉ，ｙｉ）の軌跡のうちの１つが、２次元座標のＸ軸のプラス側とマイナス側に跨っており、第１変数ｘｉの値がうまく分岐できていないことが分かる。したがって、ユーザはこのモニタリング画像を確認することで、パラメータの設定が適切でないことを確認できる。
【００７７】
なお、以上で説明した例では、モニタリング画像を静止画として表示装置１１０に表示させることを想定し、予め定めたステップ数分の第１変数ｘｉおよび第２変数ｙｉの値が蓄積されるたびにモニタリング画像を更新するものとしたが、これに限らない。例えば、処理部２０が第１変数ｘｉおよび第２変数ｙｉの値を出力するたびに座標点（ｘｉ，ｙｉ）を２次元座標上にプロットし、座標点（ｘｉ，ｙｉ）の軌跡をアニメーションとして表すモニタリング画像を表示装置１１０に表示させるようにしてもよい。この場合、軌跡を構成する座標点（ｘｉ，ｙｉ）の数に上限を設け、軌跡を構成する座標点（ｘｉ，ｙｉ）の数が上限に達したら、新たな座標点（ｘｉ，ｙｉ）が追加される際に最も古い座標点（ｘｉ，ｙｉ）を消去することで、座標点（ｘｉ，ｙｉ）の軌跡を見やすくすることが望ましい。
【００７８】
また、以上で説明した例では、処理部２０が上述の処理手順（ループ処理）を繰り返している間に、表示制御部３０がモニタリング画像を生成して表示装置１１０に表示させることを想定したが、これに限らない。例えば、表示制御部３０は、処理部２０が上述の処理手順を繰り返している間、１ステップごとに出力される第１変数ｘｉおよび第２変数ｙｉの値を蓄積しておき、上述の処理手順の繰り返しが終了した後に、蓄積した第１変数ｘｉおよび第２変数ｙｉの値に基づいて、上述のような座標点（ｘｉ，ｙｉ）の軌跡を静止画あるいはアニメーションとして表すモニタリング画像を表示装置１１０に表示させるようにしてもよい。
【００７９】
以上のように第１変数ｘｉおよび第２変数ｙｉの時間変化を座標点（ｘｉ，ｙｉ）の軌跡として表すモニタリング画像は、「Ｎ」の値、つまりスピン数が多くなると多数の軌跡が重ね合されて、視認性が低下する。そこで、表示制御部３０は、第１変数ｘｉおよび第２変数ｙｉの時間変化を濃度で表すモニタリング画像を生成し、表示装置１１０に表示させるようにしてもよい。
【００８０】
例えば、表示制御部３０は、上述のループ処理の１ステップごとに処理部２０が出力する第１変数ｘｉおよび第２変数ｙｉの値により表される座標点（ｘｉ，ｙｉ）を、複数の小ブロックにブロック分割された２次元座標上にプロット（配置）する。そして、例えば図１５および図１６に示すように、各小ブロックをブロック内プロット数に応じた濃度で表したモニタリング画像を生成し、表示装置１１０に表示させる。なお、第１変数ｘｉおよび第２変数ｙｉは実数値のため、座標点（ｘｉ，ｙｉ）の位置がモニタリング画像として表示する２次元座標の範囲よりも外に出ることもある。この場合は、２次元座標の範囲外に出た座標点（ｘｉ，ｙｉ）を、例えばＸ軸方向の最端部に位置する小ブロックにプロットすればよい。
【００８１】
図１５および図１６に示すモニタリング画像は、「Ｎ」の値（スピン数）が１０００の場合の例であり、それぞれ、（ａ）が上述の新規な運動方程式の時間発展の初期状態に対応し、（ｂ）は時間発展がある程度進んだ状態に対応し、（ｃ）は時間発展がさらに進んだ状態に対応する。
【００８２】
図１５に示すモニタリング画像では、上述の新規な運動方程式の時間発展が進むにつれて、２次元座標のＸ軸のプラス側とマイナス側の一定値付近にある小ブロックの濃度が高くなっており、第１変数ｘｉの値がうまく分岐していることが分かる。したがって、ユーザはこのモニタリング画像を参照することで、パラメータの設定が適切であることを確認できる。
【００８３】
一方、図１６に示すモニタリング画像では、上述の新規な運動方程式の時間発展が進むにつれて、２次元座標のＸ軸の両端部にある小ブロックの濃度が高くなっており、第１変数ｘｉの値が発散してしまっていることが分かる。したがって、ユーザはこのモニタリング画像を確認することで、パラメータの設定が適切でないことを確認できるとともに、第１変数ｘｉの値が発散しないように時間ステップ「ｄｔ」を小さくする必要があるといったように、パラメータのチューニングの方向性を判断することができる。
【００８４】
なお、以上で説明した例では、処理部２０が出力する第１変数ｘｉおよび第２変数ｙｉの値により表される座標点（ｘｉ，ｙｉ）を２次元座標上にプロットすることによりモニタリング画像を生成したが、これに限らない。モニタリング画像は、処理部２０が上述の新規な運動方程式の時間発展を計算するために処理手順を繰り返している間の第１変数ｘｉおよび第２変数ｙｉの時間変化を表す画像であればよく、上述の例とは異なる形態の画像であってもよい。
【００８５】
例えば、表示制御部３０は、図１７に示すように、処理部２０が出力する第１変数ｘｉの少なくとも一部の値の時間変化を表す（一つの軸が時刻である）第１グラフ（図１７の上段のグラフ）と、処理部２０が出力する第２変数ｙｉの少なくとも一部の値の時間変化を表す（一つの軸が時刻である）第２グラフ（図１７の下段のグラフ）と、の少なくともいずれかを含むモニタリング画像を表示装置１１０に表示させるようにしてもよい。図１７に示す２つのグラフは横軸が時間軸であり、第１グラフの縦軸が第１変数ｘｉの値、第２グラフの縦軸は第２変数ｙｉの値をそれぞれ示している。図１７に示す例では、時間の経過に伴って第１変数ｘｉの値がプラス側とマイナス側の一定値付近に収束していく様子が分かる。したがって、ユーザはこのモニタリング画像を参照することで、パラメータの設定が適切であることを確認できる。
【００８６】
以上、具体的な例を挙げながら詳細に説明したように、本実施形態に係る計算装置１００は、上述の新規な運動方程式の時間発展を計算することにより、イジング問題を解く構成である。新規な運動方程式は、例えばシンプレクティック・オイラー法などの計算量が小さく安定な離散解法を適用し、並列計算により高速に解くことができる。したがって、本実施形態に係る計算装置１００によれば、大規模なイジング問題を高速に解くことができる。
【００８７】
また、本実施形態に係る計算装置１００は、上述の新規な運動方程式の時間発展を計算している間の第１変数ｘｉおよび第２変数ｙｉの時間変化を表すモニタリング画像を表示装置１１０に表示させるようにしている。したがって、ユーザは、表示装置１１０に表示されているモニタリング画像を参照することにより、第１変数ｘｉの分岐が適切に行われているか否かを確認することができ、パラメータの設定が適切であるかが確認できる。また、分岐が適切でない場合に、モニタリング画像で表される第１変数ｘｉおよび第２変数ｙｉの時間変化の様子から、パラメータのチューニングの方向性を判断することができる。
【００８８】
本実施形態に係る計算装置１００の機能的な各部（上述の設定部１０、処理部２０、表示制御部３０）は、一例として、プログラム（ソフトウェア）をコンピュータに実行させることにより実現できる。コンピュータに上述の設定部１０、処理部２０、表示制御部３０などの機能を実現させるためのプログラムは、例えば、磁気ディスク（フレキシブルディスク、ハードディスクなど）、光ディスク（ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ−ＲＯＭ、ＤＶＤ±Ｒ、ＤＶＤ±ＲＷ、Ｂｌｕ−ｒａｙ（登録商標）Ｄｉｓｃなど）、半導体メモリ、またはこれに類する記録媒体に記録されて提供される。なお、プログラムを記録する記録媒体は、コンピュータが読み取り可能な記録媒体であれば、その記憶形式は何れの形態であってもよい。また、上記プログラムを、コンピュータに予めインストールするように構成してもよいし、ネットワークを介して配布される上記のプログラムをコンピュータに適宜インストールするように構成してもよい。
【００８９】
コンピュータで実行されるプログラムは、例えば、本実施形態に係る計算装置１００の機能的な各部（設定部１０、処理部２０、表示制御部３０）を含むモジュール構成となっており、実際のハードウェアとしては、例えば、第１プロセッサ回路１０１や第２プロセッサ回路１０２が上記プログラムを読み出して実行することにより、上記各部が主記憶装置１０３上にロードされ、主記憶装置１０３上に生成されるようになっている。
【００９０】
以上、本発明の実施形態を説明したが、この実施形態は例として提示したものであり、発明の範囲を限定することは意図していない。この新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
【符号の説明】
【００９１】
１０  設定部
２０  処理部
３０  表示制御部
１００  計算装置
１０１  第１プロセッサ回路
１０２  第２プロセッサ回路
１１０  表示装置
ｘｉ  第１変数
ｙｉ  第２変数
【図１】
2020154524.tif 000012
【図２】
2020154524.tif 000013
【図３】
2020154524.tif 000014
【図４】
2020154524.tif 000015
【図５】
2020154524.tif 000016
【図６】
2020154524.tif 000017
【図７】
2020154524.tif 000018
【図８】
2020154524.tif 000019
【図９】
2020154524.tif 000020
【図１０】
2020154524.tif 000021
【図１１】
2020154524.tif 000022
【図１２】
2020154524.tif 000023
【図１３】
2020154524.tif 000024
【図１４】
2020154524.tif 000025
【図１５】
2020154524.tif 000026
【図１６】
2020154524.tif 000027
【図１７】
2020154524.tif 000028