JP7228425B2 - 計算装置、表示装置およびプログラム - Google Patents

Description

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

組み合わせ最適化問題の多くは、統計力学のイジングモデルのエネルギー最小状態（基底状態）を求める問題（以下、「イジング問題」と呼ぶ）に変換できる。近年、イジング問題を高速に解く量子コンピュータ（以下、「イジングマシン」と呼ぶ）の開発が盛んに進められている。しかし、現在提案されているイジングマシンはスピン数（量子ビットの数）に限界があり、大規模なイジング問題を解くことが難しいといった課題がある。

一方、計算機によるシミュレーションでイジング問題を解く方法として、シミュレーテッド・アニーリングが知られている。シミュレーテッド・アニーリングは、大規模なイジング問題を扱うことが可能であるが、逐次更新アルゴリズムを採用するため並列計算に不向きであり、大規模なイジング問題を高速に解くことが難しい。

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

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_ｉ，ｊ」（第１パラメータ）を含む。また、パラメータは、外場を表すベクトル「ｈ_ｉ」（第２パラメータ）をさらに含んでもよい。

次に、設定部１０は、第１変数群｛ｘ｝に含まれる各第１変数ｘ_ｉおよび第２変数群｛ｙ｝に含まれる各第２変数ｙ_ｉを、例えば絶対値が０．１以下の乱数によってランダムに初期化する（ステップＳ２）。

次に、処理部２０が、後述の新規な運動方程式の時間発展の計算を行う（ステップＳ３）。すなわち、処理部２０は、上述の第１変数更新と第２変数更新とを含む処理手順を複数回繰り返すことにより、後述の新規な運動方程式の時間発展を計算する。第１変数更新では、第２変数群｛ｙ｝を用いて第１変数群｛ｘ｝が更新され、第２変数更新では、第１変数群｛ｘ｝を用いて第２変数群｛ｙ｝が更新される。複数回の処理手順の少なくとも一部において、第１変数更新の後に第２変数更新が行われてもよいし、第２変数更新の後に第１変数更新が行われてもよい。第１変数更新と第２変数更新とを含む処理手順は、所定の終了条件が満たされるまで繰り返される。

処理手順の繰り返しが終了すると、処理部２０は、最後に更新された第１変数群｛ｘ｝に含まれる各第１変数ｘ_ｉの符号を出力する（ステップＳ４）。これら第１変数ｘ_ｉの符号は、基底状態のイジングスピンを表している。本実施形態に係る計算装置１００は、以上のように、新規な運動方程式の時間発展を計算することにより、大規模なイジング問題を高速に解くことができる。

ここで、本実施形態に係る計算装置１００が扱うイジング問題について説明する。イジングエネルギーは、例えば下記式（１）のように表される。

この式（１）において、「Ｎ」はイジングスピンの数（スピン数）である。「ｓ_ｉ」はｉ番目のイジングスピンであり、例えば「ｓ_ｉ」＝±１である。「J_ｉ，ｊ」はイジングスピン間の相互作用の強さを表す行列（第１パラメータ）であり、例えば、対角成分（対角要素）が全てゼロの実対象行列である。「ｈ_ｉ」は外場を表すベクトル（第２パラメータ）である。

上記式（１）に関して、量子分岐マシンの古典モデル（以下、「古典分岐マシン」と呼ぶ）が提案されている。例えば特許文献１には、古典分岐マシンによるシミュレーションの一例が記載されている。古典分岐マシンにおいては、運動方程式は、以下の式（２）乃至式（４）の連立常微分方程式で与えられる。

これら式（２）乃至式（４）において、「Ｎ」は、例えば、イジングスピンの数（スピン数）に対応する。「Ｄ」は、例えば、特許文献１に記載の「離調」に対応するパラメータである。「ｃ」は定数である。「ｐ」は、例えば、特許文献１に記載の「ポンプ振幅（ポンプレート）」に対応するパラメータである。「Ｋ」は、例えば、特許文献１に記載の「カー係数」に対応するパラメータである。上記式（２）乃至式（４）で表される運動方程式において、第２パラメータ「ｈ_ｉ」は設けられなくてもよい。その場合は、上記式（３）および式（４）の中の「ｈ_ｉ」を含む項は無視される。

上記式（２）乃至式（４）において、ｐ（ｔ）をゼロから十分大きな値へ増加させたときの「ｘ_ｉ」の最終値（最後に更新された値）の符号「±１」が、最適解（基底状態）のイジングスピン「ｓ_ｉ」となる。「ａ（ｔ）」は、「ｐ（ｔ）」とともに増加するパラメータである。「ａ（ｔ）」は、例えば、下記式（５）により表される。

上記の古典分岐マシンは、上記式（２）乃至式（４）において、「Ｈ」をハミルトニアンとしたハミルトン力学系であると見なすことができる。

計算機によるシミュレーションでイジング問題を解く方法としては、シミュレーテッド・アニーリングが広く知られている。この方法では、逐次更新アルゴリズムが採用される。逐次更新アルゴリズムにおいては、複数のスピンが１つずつ更新される。このような逐次更新アルゴリズムは、並列計算に向かない。

これに対して、上記の古典分岐マシンの運動方程式をデジタル計算機で離散解法によって解くことが考えられる。このアルゴリズムは、シミュレーテッド・アニーリングとは異なり、並列更新アルゴリズムである。並列更新アルゴリズムにおいては、複数の変数を同時に更新することができる。このため、並列計算による高速化が期待できる。

上記式（２）乃至式（４）を用いるアプローチには、以下の課題があると考えられる。すなわち、最も計算量が大きい行列J_ｉ，ｊに関する積和計算が、第１変数ｘ_ｉの更新と第２変数ｙ_ｉの更新の双方において必要となる。上記の運動方程式は、数値的に容易には解けないため、例えば、計算量が大きい離散解法（例えば４次のルンゲ・クッタ法など）が必要となる。このため、並列計算を導入しても十分な高速化が図れない懸念がある。

そこで、本実施形態においては、上記式（２）乃至式（４）で表される運動方程式ではなく、例えば、下記式（６）乃至式（８）で表される新規な運動方程式を解く。

これら式（６）乃至式（８）において、「Ｎ」は、例えば、イジングスピンの数（スピン数）に対応する。「Ｄ」は、例えば、特許文献１に記載の「離調」に対応するパラメータである。「ｃ」は定数である。「ｐ」は、例えば、特許文献１に記載の「ポンプ振幅（ポンプレート）」に対応するパラメータ（演算パラメータ）である。「Ｋ」は、例えば、特許文献１に記載の「カー係数」に対応するパラメータである。これらのパラメータは、例えば、図３のステップＳ１において設定される。１つの例において、「Ｋ」は１に設定される。上記式（６）乃至式（８）で表される運動方程式において、第２パラメータ「ｈ_ｉ」は設けられなくてもよい。その場合は、上記式（７）および式（８）の中の「ｈ_ｉ」を含む項は無視される。

上記式（６）乃至式（８）で表される新規な運動方程式を用いる場合、最も計算量が大きい行列J_ｉ，ｊに関する積和演算は、第２変数ｙ_ｉの更新においてのみ行われ、第１変数ｘ_ｉの更新では行われない。したがって、計算量が削減される。この新規な運動方程式では、第１変数ｘ_ｉの時間微分から「ｐ」が消去されている。この新規な運動方程式において、第１変数ｘ_ｉの時間微分は、第２変数ｙ_ｉを含み、第１変数ｘ_ｉを含まない。また、第２変数ｙ_ｉの時間微分は、第１変数ｘ_ｉを含み、第２変数ｙ_ｉを含まない。ハミルトニアンにおいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉは、互いに分離されている。このため、計算量が小さく安定な離散解法が適用可能となる。例えば、シンプレクティック・オイラー法と呼ばれる方法が適用できる。

イジング問題を解くアプローチとして新規な運動方程式を用いた場合に、高い性能（例えば、高い精度）が維持できることが分かった。本実施形態に係る計算装置１００では、上記の分離可能なハミルトニアンを持つハミルトン力学系（新しい古典分岐マシン）の新規な運動方程式が、例えば、シンプレクティック・オイラー法で解かれる。本実施形態に係る計算装置１００は、このような新しいアルゴリズムの計算が並列計算によって、できる限り高速に実行されるように構成される。

実施形態においては、図３のステップＳ３において、処理部２０が第１変数更新および第２変数更新を含む処理手順を繰り返すことで上述の新規な運動方程式の時間発展を計算し、イジングモデルの基底状態を求める。例えば、「ｐ（ｔ）」をゼロから十分大きな値へ増加させたときの第１変数ｘ_ｉの最終値の符号「±１」が、基底状態のイジングスピンｓ_ｉとなる。以下では、図３のステップＳ３における処理の詳細について説明する。

図４は、処理部２０による処理の一例を示すフローチャートである。この図４のフローチャートで示す処理は、図３のステップＳ３においてサブルーチンとして実行される。

このサブルーチンが開始されると、処理部２０は、まず、時刻「ｔ」と演算パラメータ「ｐ」および「ａ」を初期化する（ステップＳ１０１）。１つの例において、時刻「ｔ」と演算パラメータ「ｐ」および「ａ」は、０に初期化される。なお、「Ｔ」は、時刻「ｔ」の最終値に対応する。「ｄｔ」は、時刻「ｔ」の１ステップ当たりの増加量である。「ｄｐ」は、演算パラメータ「ｐ」の１ステップ当たりの増加量である。これら「Ｔ」、「ｄｔ」、「ｄｐ」は、例えば、図３のステップＳ１において設定されるパラメータに含まれる。

処理部２０は、「ｔ」が「Ｔ」よりも小さいときに、以下に説明する第１変数更新および第２変数更新を含む処理手順を実施する（ステップＳ１０２）。すなわち、第１変数更新および第２変数更新を含む処理手順は、「ｔ」が「Ｔ」以上という終了条件を満たすまで繰り返されるループ処理の１ステップとして実施される。なお、「ｐ」が予め適切に設定された「Ｐ」以上であることを終了条件とし、「ｐ」が「Ｐ」よりも小さいときに、第１変数更新および第２変数更新を含む処理手順を実施するようにしてもよい。

ループ処理の１ステップとして実施される処理手順において、処理部２０は、例えば、まず、第１変数更新を行う（ステップＳ１１０）。第１変数更新では、Ｎ個のイジングスピンｓ_ｉに対応するＮ個の第１変数ｘ_ｉの更新が行われる。例えば、更新前の第１変数ｘ_ｉに、ｄｔ＊Ｄ＊ｙ_ｉを加えて得られる値を、更新後の第１変数ｘ_ｉとする。ここで、「＊」は、積の記号である。

次に、処理部２０は、第２変数更新を行う（ステップＳ１２０）。第２変数更新では、Ｎ個の第１変数ｘ_ｉに対応するＮ個の第２変数ｙ_ｉの更新が行われる。第２変数更新は、第１サブ更新（ステップＳ１２１）と第２サブ更新（ステップＳ１２２）とを含む。第１サブ更新と第２サブ更新の順番は、入れ替えが可能である。第１サブ更新の少なくとも一部と第２サブ更新の少なくとも一部とが同時に行われてもよい。

第１サブ更新では、例えば、更新前の第２変数ｙ_ｉに、ｄｔ＊［（ｐ－Ｄ－ｘ_ｉ＊ｘ_ｉ）＊ｘ_ｉ－ｃ＊ｈ_ｉ＊ａ］を加えて得られる値を、更新後の第２変数ｙ_ｉとする。

第２サブ更新では、例えば、更新前の第２変数ｙ_ｉに、ｄｔ＊ｃ＊Σ（Ｊ_ｉ，ｊ＊ｘ_ｊ）を加えて得られる値を更新後の第２変数ｙ_ｉとする。「Σ」は、行列Ｊに関する和を表す。例えば、「ｄｔ＊ｃ＊Ｊ」を行列Ｊとしてもよい。この場合、「ｄｔ＊ｃ＊」の演算を実際に行わなくてもよい。

第１変数更新および第２変数更新が終了すると、処理部２０は、更新後の第１変数ｘ_ｉおよび第２変数ｙ_ｉの値を出力する（ステップＳ１３０）。これら更新後の第１変数ｘ_ｉおよび第２変数ｙ_ｉの値は、表示制御部３０が後述のモニタリング画像を生成して表示装置１１０に表示させるために利用される。

次に、処理部２０は、時刻「ｔ」と演算パラメータ「ｐ」および「ａ」を更新する（ステップＳ１４０）。例えば、更新前の「ｔ」に「ｄｔ」を加えて得られる値を、更新後の「ｔ」とする。また、更新前の「ｐ」に「ｄｐ」を加えて得られる値を、更新後の「ｐ」とする。また、「ｐ^１／２」を、更新後の「ａ」とする。本実施形態では、「ｐ」の更新において、任意の増加関数を用いた線形増加が適用される。

そして、更新後の「ｔ」が「Ｔ」よりも小さければ、ステップＳ１０２に戻って以降の処理手順を繰り返す（ステップＳ１０３）。その後、更新後の「ｔ」が「Ｔ」以上になると、ループ処理を終了して図３のステップＳ４に進む。なお、更新後の「ｐ」が予め適切に設定された「Ｐ」よりも小さければステップＳ１０２に戻り、更新後の「ｐ」が「Ｐ」以上になると、ループ処理を終了して図３のステップＳ４に進むようにしてもよい。

なお、以上の説明では、ループ処理の各ステップにおいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉの更新が、第１変数更新（ステップＳ１１０）、第２変数更新の第１サブ更新（ステップＳ１２１）、第２変数更新の第２サブ更新（ステップＳ１２２）の順番で行われるものとしたが、これらの処理の順番は任意に変更することができる。例えば、ループ処理の少なくとも一部のステップにおいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉの更新を、図５に示すように、第２変数更新の第１サブ更新（ステップＳ１２１）、第２変数更新の第２サブ更新（ステップＳ１２２）、第１変数更新（ステップＳ１１０）の順番で行うようにしてもよい。

また、ループ処理の少なくとも一部のステップにおいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉの更新を、図６に示すように、第２変数更新の第１サブ更新（ステップＳ１２１）、第１変数更新（ステップＳ１１０）、第２変数更新の第２サブ更新（ステップＳ１２２）の順番で行うようにしてもよい。また、ループ処理の少なくとも一部のステップにおいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉの更新を、図７に示すように、第２変数更新の第２サブ更新（ステップＳ１２２）、第１変数更新（ステップＳ１１０）、第２変数更新の第１サブ更新（ステップＳ１２１）の順番で行うようにしてもよい。

本実施形態に係る計算装置１００で実行されるアルゴリズムは、例えば、以下を含む。すなわち、行列Ｊ_ｉ，ｊを取得、または、計算により定める。このとき、ベクトルｈ_ｉをさらに取得、または、計算により定めるようにしてもよい。行列Ｊ_ｉ，ｊおよびベクトルｈ_ｉは、イジングモデルにおけるパラメータである。２種類の変数（第１変数群｛ｘ｝および第２変数群｛ｙ｝）を用いる。一方の変数の更新には、他方の変数の値を用いる。一方の変数の更新において、その一方の変数の値を用いない。一方の変数の更新の後に、更新後のその一方の変数の値を用いて、他方を更新する。

すなわち、第１変数群｛ｘ｝を構成する各第１変数ｘ_ｉの更新（第１変数更新）は、更新前の第１変数ｘ_ｉに第１関数より求まる値を加えて第１変数ｘ_ｉを更新することを含む。また、第２変数群｛ｙ｝を構成する各第２変数ｙ_ｉの更新（第２変数更新）は、更新前の第２変数ｙ_ｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙ_ｉを更新することを含む。第１関数は、第１変数群｛ｘ｝から独立している。第１変数群｛ｘ｝の値を変更しても、第１関数の結果の値は変化しない。第２関数は、第２変数群｛ｙ｝から独立している。第２変数群｛ｙ｝の値を変更しても、第２関数の結果の値は変化しない。第３関数は、第２変数群｛ｙ｝から独立している。第２変数群｛ｙ｝の値を変更しても、第３関数の結果の値は変化しない。

第１関数は、例えば、ｄｔ＊Ｄ＊ｙ_ｉである（図４参照）。第２関数は、例えば、ｄｔ＊［（ｐ－Ｄ－ｘ_ｉ＊ｘ_ｉ）＊ｘ_ｉ－ｃ＊ｈ_ｉ＊ａ］である（図４参照）。第３関数は、例えば、ｄｔ＊ｃ＊Σ（Ｊ_ｉ，ｊ＊ｘ_ｊ）である（図４参照）。第１関数の変数は、第２変数ｙ_ｉを含む。第２関数の変数は、第１変数ｘ_ｉを含む。第３関数は、行列Ｊ_ｉ，ｊと第１変数ｘ_ｊとの積和演算Σ（Ｊ_ｉ，ｊ＊ｘ_ｊ）を含む。

第２関数は、例えば、第１変数ｘ_ｉの非線形関数である第４関数を含む。第４関数は、演算パラメータ「ｐ」も含む。２種類の変数の更新とともに「ｐ」は変化する。２種類の変数の更新とともに「ｐ」が変化すると、第４関数の実数根の数が変化する。「関数の実数根」とは、関数の値がゼロとなる変数の値（実数）のことである。第２変数更新において第４関数のみを考慮した場合、第４関数の実数根は非線形力学系の固定点（fixed point）に対応する。（ハミルトン力学系では、固定点はハミルトニアンの極値に対応する。）よって、第４関数の実数根の数が変化することは、固定点の数が変化することに対応する。これは、非線形力学系の分岐現象に対応する。

本実施形態に係る計算装置１００で用いるアルゴリズムでは、第１変数ｘ_ｉおよび第２変数ｙ_ｉの初期値が、初期の１つの安定固定点（stable fixed point）の近傍に設定される。「ｐ」を変化させることで分岐（bifurcation）を起こす。分岐後の複数の安定な固定点（変数の値はこの複数の安定固定点のうちの１つの近傍へと変化する）と、解きたい組合せ最適化問題の離散変数とを対応させる。これにより、分岐現象によって、その組合せ最適化問題を解く。例えば、上記の例では、分岐後の安定固定点における各第１変数ｘ_ｉの値が、正負の２値であり、その符号と、イジングスピン（イジング問題の離散変数）とが対応付けられている。初期の安定固定点が原点であるため、各第１変数ｘ_ｉと各第２変数ｙ_ｉの初期値を原点の近傍（つまり、絶対値が０．１以下の小さな乱数）の値に設定する。

第４関数は、例えば、ｄｔ＊（ｐ－Ｄ’－ｘ_ｉ＊ｘ_ｉ）＊ｘ_ｉである。「Ｄ’」は、０≦Ｄ’≦Ｄを満たす適切な定数である。初期時刻ではｐ＝０であり、第４関数の根はｘ_ｉ＝０の１つだけであるが、ｐがＤ’よりも大きくなると、根は３つとなり、そのうちの正負２つの根が、イジングスピンと対応付けられる。なお、第２関数ｄｔ＊［（ｐ－Ｄ－ｘ_ｉ＊ｘ_ｉ）＊ｘ_ｉ－ｃ＊ｈ_ｉ＊ａ］は、ｄｔ＊（ｐ－Ｄ’－ｘ_ｉ＊ｘ_ｉ）＊ｘ_ｉ＋ｄｔ＊［－（Ｄ－Ｄ’）＊ｘ_ｉ－ｃ＊ｈ_ｉ＊ａ］のように、第４関数と１次関数との和で表すことができる。したがって、第２関数は第４関数を含む。

上記の第４関数は、例えば、３次関数である。このような処理により、例えば、ニューラルネットワークで用いられている非線形関数（例えばシグモイド関数）を用いた計算よりも、計算が容易になる。

本実施形態に係る計算装置１００においては、時間ステップ（例えば「ｄｔ」）を大きくすると、計算が高速になる。一方、時間ステップを過度に大きくすると、計算が不安定になる。このことを考慮して、計算の一部における時間ステップを大きくし、計算の別の一部では、時間ステップを小さくてもよい。例えば、計算量の大きい、行列Ｊ_ｉ，ｊと第１変数ｘ_ｊとの積和演算を含む第２変数更新の第２サブ更新（ステップＳ１２２）には大きい時間ステップを適用し、その他の第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）には小さい時間ステップを適用するようにしてもよい。これにより、さらなる高速化が可能である。

このような計算を行う場合の処理部２０による処理の一例を図８に示す。図８に示す例では、１つのループ（ステップＳ１０２～ステップＳ１０３）の中に、小さいループ（ステップＳ１０４～ステップＳ１０５）が設けられている。ステップＳ１０４において、ループ変数「ｍ」は、１以上Ｍ以下である。小さいループの中において、第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）がＭ回繰り返される。これら第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）をＭ回繰り返した後、ステップＳ１２２に進んで第２変数更新の第２サブ更新が行われる。

なお、小さいループ中における第１変数更新（ステップＳ１１０）と第２変数更新の第１サブ更新（ステップＳ１２１）の順序は入れ替えが可能である。第１変数更新（ステップＳ１１０）と第２変数更新の第１サブ更新（ステップＳ１２１）の順序を入れ替えた例を図９に示す。

なお、図８および図９の例では、第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）をＭ回繰り返した後に第２変数更新の第２サブ更新（ステップＳ１２２）を行うようにしているが、図１０および図１１に示すように、第２変数更新の第２サブ更新（ステップＳ１２２）を行った後に、第１変数更新（ステップＳ１１０）および第２変数更新の第１サブ更新（ステップＳ１２１）をＭ回繰り返すようにしてもよい。

本実施形態に係る計算装置１００おいて、処理部２０は、主記憶装置１０３などのメモリに保持されたデータを読み出し、データを更新し、更新されたデータをメモリに保持させる。例えば、上述の第１変数更新（ステップＳ１１０）において、処理部２０は、更新前の第１変数ｘ_ｉと第２変数ｙ_ｉをメモリから取得し、第２変数ｙ_ｉを用いて第１関数を計算し、得られた値を更新前の第１変数ｘ_ｉに加えることにより第１変数ｘ_ｉを更新し、更新された第１変数ｘ_ｉをメモリに保持させる。

また、上述の第２変数更新の第１サブ更新（ステップＳ１２１）において、処理部２０は、更新前の第２変数ｙ_ｉと第１変数ｘ_ｉをメモリから取得し、第１変数ｘ_ｉを用いて第２関数を計算し、得られた値を更新前の第２変数ｙ_ｉに加えることにより第２変数ｙ_ｉを更新し、更新された第２変数ｙ_ｉをメモリに保持させる。また、上述の第２変数更新の第２サブ更新（ステップＳ１２２）において、処理部２０は、更新前の第２変数ｙ_ｉと行列Ｊ_ｉ，ｊと第１変数ｘ_ｊをメモリから取得し、行列Ｊ_ｉ，ｊと第１変数ｘ_ｊを用いて第３関数を計算し、得られた値を更新前の第２変数ｙ_ｉに加えることにより第２変数ｙ_ｉを更新し、更新された第２変数ｙ_ｉをメモリに保持させる。

実施形態において、例えば、行列Ｊ_ｉ，ｊが疎行列の場合、疎行列の圧縮形式を用いてもよい。疎行列の圧縮形式には、例えば、ＣＯＯ（coordinate）形式またはＣＳＲ（compressed sparse row）形式などが適用できる。疎行列の圧縮形式を用いることで、例えば、メモリサイズを節約できる。疎行列の圧縮形式を用いることで、例えば、行列Ｊ_ｉ，ｊと第１変数ｘ_ｊとの積和演算を高速に実施できる。

以下、定数「ｃ」の例について述べる。例えば、離調「Ｄ」は、行列Ｊ_ｉ，ｊの最大固有値λｍａｘのｃ倍よりも大きくされる（例えば、非特許文献２参照）。「Ｄ」が大きすぎると、無駄な計算時間が生じる。このため、例えば、「Ｄ」がλｍａｘのｃ倍と実質的に等しいように設定する。この場合、ｃ＝Ｄ／λｍａｘとなる。一方、１つの例において、行列Ｊ_ｉ，ｊは、実対称行列である。この場合に、行列Ｊ_ｉ，ｊのサイズが十分大きい場合、λｍａｘは、２σ×Ｎ１／２と実質的に同じになる。この関係は、ランダム行列のウィグナーの半円則に基づく。「σ」は、行列Ｊの非対角成分の標準偏差である。この場合、ｃ＝Ｄ／（２σ×Ｎ１／２）と設定するとよい。

本実施形態において、精度を向上させる方法として、上記式（７）の非線形関数として用いられる関数を変更することが考えられる。例えば、上記式（７）の代わりに下記式（９）を用いてもよい。

この式（９）において、「ｎ」は、２以上の偶数である。このような関数を用いることで、例えば、イジング問題の解の精度を向上することができる。

次に、本実施形態に係る計算装置１００によるモニタリング画像の表示制御について説明する。モニタリング画像は、処理部２０が上述の新規な運動方程式の時間発展を計算するために処理手順を繰り返している間の第１変数ｘ_ｉおよび第２変数ｙ_ｉの時間変化を可視化した画像である。モニタリング画像は、上述のループ処理の１ステップごとに処理部２０が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値に基づいて、表示制御部３０により生成され、表示装置１１０に表示される。

表示制御部３０は、例えば、上述のループ処理の１ステップごとに処理部２０が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値を、予め定めたステップ数ｍ分蓄積する。そして、蓄積した各ステップにおける第１変数ｘ_ｉおよび第２変数ｙ_ｉの値により表される座標点（ｘ_ｉ，ｙ_ｉ）を２次元座標上にプロットすることによりモニタリング画像を生成し、このモニタリング画像を表示装置１１０に表示させる。この際、モニタリング画像上の座標点（ｘ_ｉ，ｙ_ｉ）をｉの値に応じてＮ種類に色分けして表示させることが望ましい。

このモニタリング画像では、２次元座標上における座標点（ｘ_ｉ，ｙ_ｉ）のステップごとの位置変化を軌跡として表すことができ、上述の新規な運動方程式の時間発展に伴う第１変数ｘ_ｉおよび第２変数ｙ_ｉの時間変化を視覚的に分かり易く表示することができる。ここで、色分けしない場合であっても、それぞれのｉについて過去に蓄積されたｍ個の点を表示して、それらのｍ個の点の集合を軌跡で描くようにする、または、それらｍ個の点を順に線分で繋いで軌跡として表す、という方法を採ることができる。なお、モニタリング画像の表示中、処理部２０が新たに出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値は随時蓄積され、予め定めたステップ数分蓄積されると新たにモニタリング画像が生成され、表示装置１１０が表示するモニタリング画像が更新される。モニタリング画像の更新単位となるステップ数は、上述の新規な運動方程式の時間発展の初期で第１変数ｘ_ｉが未分岐の状態と、上述の新規な運動方程式の時間発展が進んで第１変数ｘ_ｉが分岐していることが期待される状態とが混在しないように設定するとよい。

図１２乃至図１４は、モニタリング画像の一例を示す図であり、ｉ＝１～４の４つの座標点（ｘ_ｉ，ｙ_ｉ）を各ステップごとに２次元座標上にプロットしたモニタリング画像を示している。図１２に示すモニタリング画像は、上述の新規な運動方程式の時間発展の初期状態に対応し、４つの座標点（ｘ_ｉ，ｙ_ｉ）の軌跡が２次元座標の原点（０，０）付近に集中していることから、第１変数ｘ_ｉが未分岐の状態であることが分かる。

図１３に示すモニタリング画像は、上述の新規な運動方程式の時間発展が進んで、図１２に示すモニタリング画像から更新されたモニタリング画像の一例である。この図１３の例では、４つの座標点（ｘ_ｉ，ｙ_ｉ）の軌跡がすべて２次元座標のＸ軸のプラス側またはマイナス側の一定値付近にあり、上述の新規な運動方程式の時間発展が進んだことにより４つの第１変数ｘ_ｉの値がうまく分岐していることが分かる。したがって、ユーザはこのモニタリング画像を参照することで、パラメータの設定が適切であることを確認できる。

図１４に示すモニタリング画像は、上述の新規な運動方程式の時間発展が進んで、図１２に示すモニタリング画像から更新されたモニタリング画像の他の例である。この図１４の例では、４つの座標点（ｘ_ｉ，ｙ_ｉ）の軌跡のうちの１つが、２次元座標のＸ軸のプラス側とマイナス側に跨っており、第１変数ｘ_ｉの値がうまく分岐できていないことが分かる。したがって、ユーザはこのモニタリング画像を確認することで、パラメータの設定が適切でないことを確認できる。

なお、以上で説明した例では、モニタリング画像を静止画として表示装置１１０に表示させることを想定し、予め定めたステップ数分の第１変数ｘ_ｉおよび第２変数ｙ_ｉの値が蓄積されるたびにモニタリング画像を更新するものとしたが、これに限らない。例えば、処理部２０が第１変数ｘ_ｉおよび第２変数ｙ_ｉの値を出力するたびに座標点（ｘ_ｉ，ｙ_ｉ）を２次元座標上にプロットし、座標点（ｘ_ｉ，ｙ_ｉ）の軌跡をアニメーションとして表すモニタリング画像を表示装置１１０に表示させるようにしてもよい。この場合、軌跡を構成する座標点（ｘ_ｉ，ｙ_ｉ）の数に上限を設け、軌跡を構成する座標点（ｘ_ｉ，ｙ_ｉ）の数が上限に達したら、新たな座標点（ｘ_ｉ，ｙ_ｉ）が追加される際に最も古い座標点（ｘ_ｉ，ｙ_ｉ）を消去することで、座標点（ｘ_ｉ，ｙ_ｉ）の軌跡を見やすくすることが望ましい。

また、以上で説明した例では、処理部２０が上述の処理手順（ループ処理）を繰り返している間に、表示制御部３０がモニタリング画像を生成して表示装置１１０に表示させることを想定したが、これに限らない。例えば、表示制御部３０は、処理部２０が上述の処理手順を繰り返している間、１ステップごとに出力される第１変数ｘ_ｉおよび第２変数ｙ_ｉの値を蓄積しておき、上述の処理手順の繰り返しが終了した後に、蓄積した第１変数ｘ_ｉおよび第２変数ｙ_ｉの値に基づいて、上述のような座標点（ｘ_ｉ，ｙ_ｉ）の軌跡を静止画あるいはアニメーションとして表すモニタリング画像を表示装置１１０に表示させるようにしてもよい。

以上のように第１変数ｘ_ｉおよび第２変数ｙ_ｉの時間変化を座標点（ｘ_ｉ，ｙ_ｉ）の軌跡として表すモニタリング画像は、「Ｎ」の値、つまりスピン数が多くなると多数の軌跡が重ね合されて、視認性が低下する。そこで、表示制御部３０は、第１変数ｘ_ｉおよび第２変数ｙ_ｉの時間変化を濃度で表すモニタリング画像を生成し、表示装置１１０に表示させるようにしてもよい。

例えば、表示制御部３０は、上述のループ処理の１ステップごとに処理部２０が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値により表される座標点（ｘ_ｉ，ｙ_ｉ）を、複数の小ブロックにブロック分割された２次元座標上にプロット（配置）する。そして、例えば図１５および図１６に示すように、各小ブロックをブロック内プロット数に応じた濃度で表したモニタリング画像を生成し、表示装置１１０に表示させる。なお、第１変数ｘ_ｉおよび第２変数ｙ_ｉは実数値のため、座標点（ｘ_ｉ，ｙ_ｉ）の位置がモニタリング画像として表示する２次元座標の範囲よりも外に出ることもある。この場合は、２次元座標の範囲外に出た座標点（ｘ_ｉ，ｙ_ｉ）を、例えばＸ軸方向の最端部に位置する小ブロックにプロットすればよい。

図１５および図１６に示すモニタリング画像は、「Ｎ」の値（スピン数）が１０００の場合の例であり、それぞれ、（ａ）が上述の新規な運動方程式の時間発展の初期状態に対応し、（ｂ）は時間発展がある程度進んだ状態に対応し、（ｃ）は時間発展がさらに進んだ状態に対応する。

図１５に示すモニタリング画像では、上述の新規な運動方程式の時間発展が進むにつれて、２次元座標のＸ軸のプラス側とマイナス側の一定値付近にある小ブロックの濃度が高くなっており、第１変数ｘ_ｉの値がうまく分岐していることが分かる。したがって、ユーザはこのモニタリング画像を参照することで、パラメータの設定が適切であることを確認できる。

一方、図１６に示すモニタリング画像では、上述の新規な運動方程式の時間発展が進むにつれて、２次元座標のＸ軸の両端部にある小ブロックの濃度が高くなっており、第１変数ｘ_ｉの値が発散してしまっていることが分かる。したがって、ユーザはこのモニタリング画像を確認することで、パラメータの設定が適切でないことを確認できるとともに、第１変数ｘ_ｉの値が発散しないように時間ステップ「ｄｔ」を小さくする必要があるといったように、パラメータのチューニングの方向性を判断することができる。

なお、以上で説明した例では、処理部２０が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値により表される座標点（ｘ_ｉ，ｙ_ｉ）を２次元座標上にプロットすることによりモニタリング画像を生成したが、これに限らない。モニタリング画像は、処理部２０が上述の新規な運動方程式の時間発展を計算するために処理手順を繰り返している間の第１変数ｘ_ｉおよび第２変数ｙ_ｉの時間変化を表す画像であればよく、上述の例とは異なる形態の画像であってもよい。

例えば、表示制御部３０は、図１７に示すように、処理部２０が出力する第１変数ｘ_ｉの少なくとも一部の値の時間変化を表す（一つの軸が時刻である）第１グラフ（図１７の上段のグラフ）と、処理部２０が出力する第２変数ｙ_ｉの少なくとも一部の値の時間変化を表す（一つの軸が時刻である）第２グラフ（図１７の下段のグラフ）と、の少なくともいずれかを含むモニタリング画像を表示装置１１０に表示させるようにしてもよい。図１７に示す２つのグラフは横軸が時間軸であり、第１グラフの縦軸が第１変数ｘ_ｉの値、第２グラフの縦軸は第２変数ｙ_ｉの値をそれぞれ示している。図１７に示す例では、時間の経過に伴って第１変数ｘ_ｉの値がプラス側とマイナス側の一定値付近に収束していく様子が分かる。したがって、ユーザはこのモニタリング画像を参照することで、パラメータの設定が適切であることを確認できる。

以上、具体的な例を挙げながら詳細に説明したように、本実施形態に係る計算装置１００は、上述の新規な運動方程式の時間発展を計算することにより、イジング問題を解く構成である。新規な運動方程式は、例えばシンプレクティック・オイラー法などの計算量が小さく安定な離散解法を適用し、並列計算により高速に解くことができる。したがって、本実施形態に係る計算装置１００によれば、大規模なイジング問題を高速に解くことができる。

また、本実施形態に係る計算装置１００は、上述の新規な運動方程式の時間発展を計算している間の第１変数ｘ_ｉおよび第２変数ｙ_ｉの時間変化を表すモニタリング画像を表示装置１１０に表示させるようにしている。したがって、ユーザは、表示装置１１０に表示されているモニタリング画像を参照することにより、第１変数ｘ_ｉの分岐が適切に行われているか否かを確認することができ、パラメータの設定が適切であるかが確認できる。また、分岐が適切でない場合に、モニタリング画像で表される第１変数ｘ_ｉおよび第２変数ｙ_ｉの時間変化の様子から、パラメータのチューニングの方向性を判断することができる。

本実施形態に係る計算装置１００の機能的な各部（上述の設定部１０、処理部２０、表示制御部３０）は、一例として、プログラム（ソフトウェア）をコンピュータに実行させることにより実現できる。コンピュータに上述の設定部１０、処理部２０、表示制御部３０などの機能を実現させるためのプログラムは、例えば、磁気ディスク（フレキシブルディスク、ハードディスクなど）、光ディスク（ＣＤ－ＲＯＭ、ＣＤ－Ｒ、ＣＤ－ＲＷ、ＤＶＤ－ＲＯＭ、ＤＶＤ±Ｒ、ＤＶＤ±ＲＷ、Ｂｌｕ－ｒａｙ（登録商標）Ｄｉｓｃなど）、半導体メモリ、またはこれに類する記録媒体に記録されて提供される。なお、プログラムを記録する記録媒体は、コンピュータが読み取り可能な記録媒体であれば、その記憶形式は何れの形態であってもよい。また、上記プログラムを、コンピュータに予めインストールするように構成してもよいし、ネットワークを介して配布される上記のプログラムをコンピュータに適宜インストールするように構成してもよい。

コンピュータで実行されるプログラムは、例えば、本実施形態に係る計算装置１００の機能的な各部（設定部１０、処理部２０、表示制御部３０）を含むモジュール構成となっており、実際のハードウェアとしては、例えば、第１プロセッサ回路１０１や第２プロセッサ回路１０２が上記プログラムを読み出して実行することにより、上記各部が主記憶装置１０３上にロードされ、主記憶装置１０３上に生成されるようになっている。

以上、本発明の実施形態を説明したが、この実施形態は例として提示したものであり、発明の範囲を限定することは意図していない。この新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。

１０ 設定部
２０ 処理部
３０ 表示制御部
１００ 計算装置
１０１ 第１プロセッサ回路
１０２ 第２プロセッサ回路
１１０ 表示装置
ｘ_ｉ 第１変数
ｙ_ｉ 第２変数

Claims (8)

1. 第１変数更新と第２変数更新とを含む処理手順を繰り返す処理部と、
   表示装置に画像を表示させる表示制御部と、を備え、
   前記第１変数更新は、Ｎ個（Ｎは２以上の整数）の第１変数ｘ_ｉ（ｉは１以上Ｎ以下の整数）の各々について、更新前の第１変数ｘ_ｉに第１関数より求まる値を加えて第１変数ｘ_ｉを更新することを含み、
   前記第２変数更新は、Ｎ個の第１変数ｘ_ｉに対応するＮ個の第２変数ｙ_ｉの各々について、更新前の第２変数ｙ_ｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙ_ｉを更新することを含み、
   前記第１関数の変数は、第２変数ｙ_ｉを含み、
   前記第２関数の変数は、第１変数ｘ_ｉを含み、
   前記第３関数は、行列J_ｉ，ｊと第１変数ｘ_ｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含み、
   前記処理部は、前記処理手順を繰り返している間は更新された第１変数ｘ_ｉおよび第２変数ｙ_ｉの値を出力し、前記処理手順の繰り返しが終了すると、最後に更新された第１変数ｘ_ｉの符号を出力し、
   前記表示制御部は、前記処理部が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値に基づいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉの少なくとも一部の時間変化を表すモニタリング画像を前記表示装置に表示させること
   を特徴とする計算装置。
2. 前記表示制御部は、前記処理部が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値により表される座標点（ｘ_ｉ，ｙ_ｉ）の少なくとも一部を２次元座標上にプロットした前記モニタリング画像を前記表示装置に表示させること
   を特徴とする請求項１に記載の計算装置。
3. 前記表示制御部は、前記処理部が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値により表される座標点（ｘ_ｉ，ｙ_ｉ）を、複数の小ブロックにブロック分割された２次元座標上にプロットし、各小ブロックをブロック内プロット数に応じた濃度で表した前記モニタリング画像を前記表示装置に表示させること
   を特徴とする請求項１に記載の計算装置。
4. 前記表示制御部は、前記プロットされた前記座標点の集合を軌跡として表すこと
   を特徴とする請求項２または３に記載の計算装置。
5. 前記表示制御部は、前記プロットされた前記座標点の間を線分で繋いで軌跡として表すこと
   を特徴とする請求項４に記載の計算装置。
6. 前記表示制御部は、前記処理部が出力する第１変数ｘ_ｉの少なくとも一部の値の時間変化を表す第１グラフと前記処理部が出力する第２変数ｙ_ｉの少なくとも一部の値の時間変化を表す第２グラフと、の少なくともいずれかを含む前記モニタリング画像を前記表示装置に表示させること
   を特徴とする請求項１に記載の計算装置。
7. 計算装置に接続される表示装置であって、
   前記計算装置は、
   第１変数更新と第２変数更新とを含む処理手順を繰り返し、
   前記第１変数更新は、Ｎ個（Ｎは２以上の整数）の第１変数ｘ_ｉ（ｉは１以上Ｎ以下の整数）の各々について、更新前の第１変数ｘ_ｉに第１関数より求まる値を加えて第１変数ｘ_ｉを更新することを含み、
   前記第２変数更新は、Ｎ個の第１変数ｘ_ｉに対応するＮ個の第２変数ｙ_ｉの各々について、更新前の第２変数ｙ_ｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙ_ｉを更新することを含み、
   前記第１関数の変数は、第２変数ｙ_ｉを含み、
   前記第２関数の変数は、第１変数ｘ_ｉを含み、
   前記第３関数は、行列J_ｉ，ｊと第１変数ｘ_ｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含み、
   前記処理手順を繰り返している間は更新された第１変数ｘ_ｉおよび第２変数ｙ_ｉの値を出力し、前記処理手順の繰り返しが終了すると、最後に更新された第１変数ｘ_ｉの符号を出力し、
   前記計算装置が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値に基づいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉの少なくとも一部の時間変化を表すモニタリング画像を表示すること
   を特徴とする表示装置。
8. コンピュータに、
   第１変数更新と第２変数更新とを含む処理手順を繰り返す処理部の機能と、
   表示装置に画像を表示させる表示制御部の機能と、を実現させるためのプログラムであって、
   前記第１変数更新は、Ｎ個（Ｎは２以上の整数）の第１変数ｘ_ｉ（ｉは１以上Ｎ以下の整数）の各々について、更新前の第１変数ｘ_ｉに第１関数より求まる値を加えて第１変数ｘ_ｉを更新することを含み、
   前記第２変数更新は、Ｎ個の第１変数ｘ_ｉに対応するＮ個の第２変数ｙ_ｉの各々について、更新前の第２変数ｙ_ｉに第２関数より求まる値および第３関数より求まる値を加えて第２変数ｙ_ｉを更新することを含み、
   前記第１関数の変数は、第２変数ｙ_ｉを含み、
   前記第２関数の変数は、第１変数ｘ_ｉを含み、
   前記第３関数は、行列J_ｉ，ｊと第１変数ｘ_ｊ（ｊはｉとは異なる１以上Ｎ以下の整数）との積和演算を含み、
   前記処理部は、前記処理手順を繰り返している間は更新された第１変数ｘ_ｉおよび第２変数ｙ_ｉの値を出力し、前記処理手順の繰り返しが終了すると、最後に更新された第１変数ｘ_ｉの符号を出力し、
   前記表示制御部は、前記処理部が出力する第１変数ｘ_ｉおよび第２変数ｙ_ｉの値に基づいて、第１変数ｘ_ｉおよび第２変数ｙ_ｉの少なくとも一部の時間変化を表すモニタリング画像を前記表示装置に表示させること
   を特徴とするプログラム。

Images