JP2019008502A - 結合振動子系の計算装置、プログラム及び方法 - Google Patents

Abstract

【課題】組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算するにあたり、光学回路を安定化する必要をなくし、機械学習のボルツマンサンプリングを行なうために、機械学習を連続値スピンモデルへマッピングし、連続値スピンモデルの平衡状態を計算するにあたり、光学回路を安定化する必要をなくす。【解決手段】結合振動子系のダイナミクスを記述する微分方程式を、光学回路で実装するのではなく、電子回路で実装し、上記の微分方程式のイテレーションの反復後に、上記の微分方程式の定常状態における各振動子毎の位相を計算することにより、結合振動子系の相互結合に対応するイジングモデル又は連続値スピンモデルの基底状態又は平衡状態を計算する。【選択図】図３

Description

本開示は、スピン系の基底状態又は平衡状態を計算することにより、そのスピン系へマッピングされた組み合わせ最適化問題の最適解を求めて、そのスピン系へマッピングされた機械学習のボルツマンサンプリングを行なう技術に関する。

組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算する技術が知られている。機械学習のボルツマンサンプリングを行なうために、機械学習を連続値スピンモデルへマッピングし、連続値スピンモデルの平衡状態を計算する技術も知られている。

特許文献１、２では、複数のレーザパルスで複数のイジングスピンを実装する。ここで、複数のレーザパルスの間の相互注入を行ない、複数のイジングスピンの間の相互結合を実装する。そして、複数のレーザパルス全体の閾値利得が最低となる発振モードを立ち上げ、複数のイジングスピン全体のエネルギーが最低となる基底状態を実装する。

特表２０１６−５２８６１１号公報 国際公開第２０１５／１５６１２６号

特許文献１のイジングモデルの計算装置Ｃ１の構成を図１に示す。レーザ発振部１１は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・を生成する。リング共振器１２は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・を周回させる。

パルス遅延線１４は、１パルス間隔分の遅延量を有する。相互結合実装部１３は、パルス遅延線１４を伝搬して相互結合係数Ｊ_１２、Ｊ_２３、Ｊ_３４、Ｊ_４１を反映した振幅及び位相を有するレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４を、パルス遅延線１４を伝搬せずリング共振器１２を伝搬するレーザパルスＰ２、Ｐ３、Ｐ４、Ｐ１にそれぞれ相互注入する。

パルス遅延線１６は、２パルス間隔分の遅延量を有する。相互結合実装部１５は、パルス遅延線１６を伝搬して相互結合係数Ｊ_１３、Ｊ_２４、Ｊ_３１、Ｊ_４２を反映した振幅及び位相を有するレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４を、パルス遅延線１６を伝搬せずリング共振器１２を伝搬するレーザパルスＰ３、Ｐ４、Ｐ１、Ｐ２にそれぞれ相互注入する。

パルス遅延線１８は、３パルス間隔分の遅延量を有する。相互結合実装部１７は、パルス遅延線１８を伝搬して相互結合係数Ｊ_１４、Ｊ_２１、Ｊ_３２、Ｊ_４３を反映した振幅及び位相を有するレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４を、パルス遅延線１８を伝搬せずリング共振器１２を伝搬するレーザパルスＰ４、Ｐ１、Ｐ２、Ｐ３にそれぞれ相互注入する。

スピン測定部１９は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・が定常状態に到達した後に、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・の位相を測定することにより、複数のレーザパルスの位相に対応する複数のイジングスピンの方向を測定する。

しかし、特許文献１では、リング共振器１２の光路長を安定化しなければ、複数のレーザパルスを安定化することができない。そして、パルス遅延線１４、１６、１８の光路長を安定化しなければ、複数のレーザパルスの間の相互注入を正確に実行することができず、複数のイジングスピンの間の相互結合を正確に実装することができない。

特許文献２のイジングモデルの計算装置Ｃ２の構成を図２に示す。レーザ発振部２１は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・を生成する。リング共振器２２は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・を周回させる。

スピン測定部２３は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・の位相を測定する。相互結合実装部２４は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・の間の相互注入を行ない、複数のイジングスピンの間の相互結合を実装する。

相互結合実装部２４は、相互結合係数ΣＪ_１ｋσ_ｋ（ｋ＝２、３、４）を電子回路上で計算し、相互結合係数ΣＪ_１ｋσ_ｋ（ｋ＝２、３、４）を反映した振幅及び位相を有するレーザパルスを生成し、生成したレーザパルスをレーザパルスＰ１に注入する。

相互結合実装部２４は、相互結合係数ΣＪ_２ｋσ_ｋ（ｋ＝１、３、４）を電子回路上で計算し、相互結合係数ΣＪ_２ｋσ_ｋ（ｋ＝１、３、４）を反映した振幅及び位相を有するレーザパルスを生成し、生成したレーザパルスをレーザパルスＰ２に注入する。

相互結合実装部２４は、相互結合係数ΣＪ_３ｋσ_ｋ（ｋ＝１、２、４）を電子回路上で計算し、相互結合係数ΣＪ_３ｋσ_ｋ（ｋ＝１、２、４）を反映した振幅及び位相を有するレーザパルスを生成し、生成したレーザパルスをレーザパルスＰ３に注入する。

相互結合実装部２４は、相互結合係数ΣＪ_４ｋσ_ｋ（ｋ＝１、２、３）を電子回路上で計算し、相互結合係数ΣＪ_４ｋσ_ｋ（ｋ＝１、２、３）を反映した振幅及び位相を有するレーザパルスを生成し、生成したレーザパルスをレーザパルスＰ４に注入する。

スピン測定部２３は、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・が定常状態に到達した後に、複数のレーザパルスＰ１、Ｐ２、Ｐ３、Ｐ４、・・・の位相を測定することにより、複数のレーザパルスの位相に対応する複数のイジングスピンの方向を測定する。

つまり、特許文献２では、相互結合実装部２４が有する電子回路を用いて、複数のレーザパルスの間の相互注入を正確に実行することができ、複数のイジングスピンの間の相互結合を正確に実装することができる。しかし、リング共振器２２の光路長を安定化しなければ、複数のレーザパルスを安定化することができない。

そこで、前記課題を解決するために、本開示は、組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算するにあたり、光学回路を安定化する必要をなくすことを目的とする。

そして、前記目的に加えて、本開示は、機械学習のボルツマンサンプリングを行なうために、機械学習を連続値スピンモデルへマッピングし、連続値スピンモデルの平衡状態を計算するにあたり、光学回路を安定化する必要をなくすことも目的とする。

前記課題を解決するために、結合振動子系のダイナミクスを記述する微分方程式を、光学回路で実装するのではなく、電子回路で実装することとした。そのうえで、上記の微分方程式のイテレーションの反復後に、上記の微分方程式の定常状態における各振動子毎の位相を計算することにより、結合振動子系の相互結合に対応するイジングモデル又は連続値スピンモデルの基底状態又は平衡状態を計算することとした。

具体的には、本開示は、結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーション部と、前記イテレーション部により反復されたイテレーションの後に、前記微分方程式の定常状態における各振動子毎の位相を計算することにより、前記結合振動子系の相互結合に対応するスピン系の基底状態又は平衡状態を計算するスピン系計算部と、を備え、前記イテレーション部は、前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算部と、前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算部と、前記微分方程式の各イテレーションにおいて、前記利得損失計算部及び前記相互結合計算部により各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算部と、を備えることを特徴とする結合振動子系の計算装置である。

また、本開示は、結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーションステップと、前記イテレーションステップにより反復されたイテレーションの後に、前記微分方程式の定常状態における各振動子毎の位相を計算することにより、前記結合振動子系の相互結合に対応するスピン系の基底状態又は平衡状態を計算するスピン系計算ステップと、をコンピュータに実行させ、前記イテレーションステップは、前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算ステップと、前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算ステップと、前記微分方程式の各イテレーションにおいて、前記利得損失計算ステップ及び前記相互結合計算ステップにより各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算ステップと、を備えることを特徴とする結合振動子系の計算プログラムである。

また、本開示は、結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーションステップと、前記イテレーションステップにより反復されたイテレーションの後に、前記微分方程式の定常状態における各振動子毎の位相を計算することにより、前記結合振動子系の相互結合に対応するスピン系の基底状態又は平衡状態を計算するスピン系計算ステップと、を備え、前記イテレーションステップは、前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算ステップと、前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算ステップと、前記微分方程式の各イテレーションにおいて、前記利得損失計算ステップ及び前記相互結合計算ステップにより各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算ステップと、を備えることを特徴とする結合振動子系の計算方法である。

この構成によれば、組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。

そして、機械学習のボルツマンサンプリングを行なうために、機械学習を連続値スピンモデルへマッピングし、連続値スピンモデルの平衡状態を計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。

また、本開示は、前記微分方程式の各イテレーションにおいて、各振動子毎の雑音に起因する各振動子毎の振幅及び位相の微小時間変化を計算する雑音計算部、をさらに備え、前記振幅位相計算部は、前記微分方程式の各イテレーションにおいて、前記利得損失計算部、前記相互結合計算部及び前記雑音計算部により各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算することを特徴とする結合振動子系の計算装置である。

この構成によれば、結合振動子系のダイナミクスを記述する微分方程式の定常状態として、複数の結合振動子の振幅が０である状態を回避することができる。

また、本開示は、結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーション部、を備え、前記イテレーション部は、前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算部と、前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算部と、前記微分方程式の各イテレーションにおいて、各振動子毎の雑音に起因する各振動子毎の振幅及び位相の微小時間変化を計算する雑音計算部と、前記微分方程式の各イテレーションにおいて、前記利得損失計算部、前記相互結合計算部及び前記雑音計算部により各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算部と、を備えることを特徴とする結合振動子系の計算装置である。

この構成によれば、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。

また、本開示は、前記雑音計算部は、前記微分方程式の初期状態のみにおいて、各振動子毎の雑音に起因する各振動子毎の振幅及び位相の微小時間変化を計算することを特徴とする結合振動子系の計算装置である。

この構成によれば、複数の結合振動子の振幅が０である初期状態から脱出することができた後には、雑音計算部の計算負担をなくすことができる。

また、本開示は、前記イテレーション部は、前回のイテレーションにおける各振動子毎の微小時間後の振幅及び位相が、前記微分方程式の定常状態に所定の精度で接近したときに、次回のイテレーションにおける各振動子毎の微小時間後の振幅及び位相の計算を停止することを特徴とする結合振動子系の計算装置である。

この構成によれば、複数の結合振動子の振幅が有限値である定常状態に接近することができた後には、イテレーション部の計算負担をなくすことができる。

また、本開示は、前記イテレーション部は、前記微分方程式の各イテレーションにおいて、各振動子毎の並列処理により、各振動子毎の微小時間後の振幅及び位相を計算することを特徴とする結合振動子系の計算装置である。

結合振動子系の全振動子数をＭとすると、各振動子毎の並列処理を行わないときには、２体相互結合の計算量は、Ｍ^２のオーダーに増えるが、各振動子毎の並列処理を行なうときには、２体相互結合の計算量は、ｌｏｇ_２Ｍのオーダーに減らせる。

また、本開示は、前記イテレーション部は、複数のイテレーション部を備え、前記複数のイテレーション部の各々は、前記微分方程式の各イテレーションにおいて、前記結合振動子系のうちの各々に特定のグループに属する振動子の微小時間後の振幅及び位相のデータを、前記複数のイテレーション部のその他との間で共有することを特徴とする結合振動子系の計算装置である。

この構成によれば、結合振動子系の全振動子数が膨大な数になるときでも、各々のイテレーション部を各々のプロセッサに実装することができる。

また、本開示は、前回のイテレーションにおける前記データの共有処理が、次回のイテレーションにおける計算処理を停止／中断させないように、前記複数のイテレーション部の各々は、前記結合振動子系のうちの各々に特定のグループに属する振動子を割り当てられることを特徴とする結合振動子系の計算装置である。

この構成によれば、各々のイテレーション部を各々のプロセッサに実装するときでも、各々のイテレーション部の間のデータの共有処理を円滑に実行することができる。

また、本開示は、前記振幅位相計算部は、前記微分方程式の各イテレーションにおいて、所定の精度で量子化したうえで、各振動子毎の微小時間後の振幅及び位相を計算することを特徴とする結合振動子系の計算装置である。

この構成によれば、結合振動子系の計算装置の計算負担を減らすことができるとともに、結合振動子系の計算装置のメモリ容量を減らすことができる。

また、本開示は、前記振幅位相計算部により計算される各振動子毎の微小時間後の振幅についての量子化カットオフは、前記振幅位相計算部により計算される各振動子毎の微小時間後の振幅についての量子化の前記所定の精度を満足するように設定されており、かつ、前記微分方程式の定常状態における各振動子毎の振幅と比べて大きく設定されていることを特徴とする結合振動子系の計算装置である。

この構成によれば、各振動子毎の振幅の分解能を低くし過ぎず、かつ、各振動子毎の振幅の制限を厳しくし過ぎず、結合振動子系の計算装置の計算負担を減らすことができるとともに、結合振動子系の計算装置のメモリ容量を減らすことができる。

また、本開示は、前記相互結合計算部により計算される各振動子毎の振幅の微小時間変化についての量子化ビット数は、前記振幅位相計算部により計算される各振動子毎の微小時間後の振幅についての量子化ビット数と比べて、少なくともｌｏｇ_２Ｍ（ただし、Ｍは、前記結合振動子系の全振動子数である。）だけ大きく設定されていることを特徴とする結合振動子系の計算装置である。

この構成によれば、各振動子毎の振幅の分解能と比べて、各振動子毎の相互結合の大きさの分解能を低くし過ぎず、結合振動子系の計算装置の計算負担を減らすことができるとともに、結合振動子系の計算装置のメモリ容量を減らすことができる。

このように、本開示によれば、組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算するにあたり、光学回路を安定化する必要をなくすことができる。

そして、本開示によれば、機械学習のボルツマンサンプリングを行なうために、機械学習を連続値スピンモデルへマッピングし、連続値スピンモデルの平衡状態を計算するにあたり、光学回路を安定化する必要をなくすことができる。

特許文献１のイジングモデルの計算装置の構成を示す図である。 特許文献２のイジングモデルの計算装置の構成を示す図である。 第１実施形態の結合振動子系の計算装置の構成を示す図である。 第１実施形態の利得損失計算部の処理内容を示す図である。 第１実施形態の結合振動子系の計算装置の計算結果を示す図である。 第１実施形態のイテレーション部の第１構成を示す図である。 第１実施形態のイテレーション部の第２構成を示す図である。 第１実施形態のイテレーション部の第３構成を示す図である。 第１実施形態の振幅位相計算部の量子化を示す図である。 第１実施形態の結合振動子系の計算装置の計算結果を示す図である。 第１実施形態の結合振動子系の計算装置の計算結果を示す図である。 第１実施形態の結合振動子系の計算装置の計算結果を示す図である。 第２実施形態の結合振動子系の計算装置の構成を示す図である。 第３実施形態の結合振動子系の計算装置の構成を示す図である。

添付の図面を参照して本開示の実施形態を説明する。以下に説明する実施形態は本開示の実施の例であり、本開示は以下の実施形態に制限されるものではない。

（第１実施形態：計算原理）
第１実施形態では、組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算する。または、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算する。

そこで、複数の振動子で複数のイジングスピンを実装する。ここで、複数の振動子の間の相互結合を行なうことにより、複数のイジングスピンの間の相互結合を実装する。そして、複数の振動子全体の閾値利得が最低となる発振モードを立ち上げることにより、複数のイジングスピン全体のエネルギーが最低となる基底状態を実装する。

まず、複数の振動子全体の閾値利得が最低となる発振モードを立ち上げることにより、複数のイジングスピン全体のエネルギーが最低となる基底状態を実装することができる理由について説明する。イジングモデルのハミルトニアンは、数１で表される。ここで、Ｊ_ｉｊは、スピンｉ、ｊの間の相互結合係数である。そして、σ_ｉ、σ_ｊ（＝＋１又は−１）は、それぞれ、スピンｉ、ｊの位相（０ｒａｄ又はπｒａｄ）を表わす。

各振動子ｉにおいて、Ｉ成分強度ｃ_ｉ及びＱ成分強度ｓ_ｉについて、レート方程式は、ファンデルポール方程式に対応して、数２、３のようになる。

ｔは、無次元時間であり、ｔ＝γ_Ｓτ／２である。τは、実時間である。γ_Ｓは、振動子の強度減衰率である。ｃ_ｉ及びｓ_ｉは、それぞれ、規格化後のＩ成分及びＱ成分の強度であり、ｃ_ｉ＝Ｃ_ｉ／Ａ_ｓ及びｓ_ｉ＝Ｓ_ｉ／Ａ_ｓである。Ｃ_ｉ及びＳ_ｉは、それぞれ、規格化前のＩ成分及びＱ成分の強度である。規格化因子Ａ_ｓは、ｐ（後述する規格化後のポンプレート）＝２における振動子の強度であり、Ａ_ｓ＝√（γ_Ｓγ_ｐ／２κ^２）である。γ_ｐは、ポンピングの強度減衰率である。κは、縮退パラメトリックゲインである。ｐは、規格化後のポンプレートであり、ｐ＝Ｆ_ｐ／Ｆ_ｔｈである。Ｆ_ｐは、規格化前のポンプレートである。規格化因子Ｆ_ｔｈは、閾値ポンプレートであり、Ｆ_ｔｈ＝γ_Ｓ√（γ_ｐ）／４κである。

数２の−ｃ_ｉ及び数３の−ｓ_ｉは、損失に関わる項である。数２の＋ｐｃ_ｉ及び数３の−ｐｓ_ｉは、線形利得に関わる項である。数２の−（ｃ_ｉ ^２＋ｓ_ｉ ^２）ｃ_ｉ及び数３の−（ｃ_ｉ ^２＋ｓ_ｉ ^２）ｓ_ｉは、飽和利得に関わる項である。これらの項は、２体の相互結合等による摂動項を含まない、ファンデルポール方程式を構成する。

数２、３のξ_ｉｊ（Ｊ_ｉｊに比例）が関わる項は、２体の相互結合に関わる項であり、ファンデルポール方程式に対する、２体の相互結合による摂動項である。数２、３のｄＷ（確率微分）が関わる項は、振動子強度の増幅段階等での雑音に関わる項であり、ファンデルポール方程式に対する、振動子強度の増幅段階等での雑音による摂動項である。

定常状態において、数２、３は、それぞれ、数４、５のようになる。

数２のｐ−（ｃ_ｉ ^２＋ｓ_ｉ ^２）は、各振動子ｉについての飽和利得である。ここで、定常状態において、複数の振動子全体についての飽和利得は、複数の振動子全体についての強度減衰率に等しく、Ｉ成分ｃ_ｉは、有限値であるが、Ｑ成分ｓ_ｉは、０である。よって、複数の振動子全体についての強度減衰率Γは、数６のようになる。

ここで、数６の最右辺の第１項は、数４の左辺の第３項を摂動項としたときにおける、摂動の０次の寄与を示す。そして、数６の最右辺の第２項は、数４の左辺の第３項を摂動項としたときにおける、摂動の１次の寄与を示す。さらに、σ_ｉ＝ｓｇｎ（ｃ_ｉ）〜ｓｇｎ（ｃ_ｉ ^（０））（ｃ_ｉ ^（０）は摂動の０次の寄与）を用いている。

ポンプレートｐを漸増制御するときには、複数の振動子全体として、「最初に」最小の強度減衰率Γを実現する発振位相状態｛σ_ｉ｝が選択される。つまり、複数の振動子全体として、１個の特定の発振モードが選択される。そして、発振モードの間の競合に起因して、１個の特定の発振モードは、他の発振モードを抑制する。つまり、複数の振動子全体として、数６のΓは最小化される。一方で、複数の振動子全体として、数６のＭは一定である。よって、複数の振動子全体として、数６の−Σξ_ｉｊσ_ｉσ_ｊは最小化される。つまり、数１のイジングハミルトニアンを最小化する基底状態が実現されたことになる。

（第１実施形態：計算方法）
次に、第１実施形態の結合振動子系の計算装置の構成を図３に示す。結合振動子系の計算装置Ｃ３は、イテレーション部３１及びスピン系計算部３２から構成される。結合振動子系の計算装置Ｃ３に対して、結合振動子系の計算プログラムをインストールすることにより、イテレーション部３１及びスピン系計算部３２を実装することができる。

イテレーション部３１は、結合振動子系のダイナミクスを記述する微分方程式（数２）の各イテレーションにおいて、各振動子ｉ毎の微小時間後の振幅及び位相（数２のｃ_ｉの絶対値及び符号）を計算する。よって、イテレーション部３１は、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算することができる。

スピン系計算部３２は、イテレーション部３１により反復されたイテレーションの後に、微分方程式（数２）の定常状態における各振動子ｉ毎の位相（数２のｃ_ｉの符号）を計算することにより、結合振動子系の相互結合（数２のξ_ｉｊ）に対応するイジングモデル（数１）の基底状態を計算する。よって、スピン系計算部３２は、イジングモデルの基底状態を計算することにより、組み合わせ最適化問題の最適解を求めることができる。

ここで、結合振動子系の全振動子数をＭとすると、数２の微分方程式のうち、右辺第１項の利得損失項及び右辺第３項の雑音項の計算量は、Ｍに比例するが、右辺第２項の相互結合項の計算量は、Ｍ^２に比例する。そこで、数２の微分方程式のうち、右辺第２項の相互結合項の計算時間が、右辺第１項の利得損失項及び右辺第３項の雑音項の計算時間と比べて、十分に長いときには、図３に示したように各イテレーションを行なえばよい。

つまり、ある回のイテレーションの後でのｃ_ｉを用いて計算した右辺第１項の利得損失項及び右辺第３項の雑音項と、その回のイテレーションよりも数回（例えば、１回のみ。）遅れた回のイテレーションの後でのｃ_ｉを用いて計算した右辺第２項の相互結合項と、を加算することにより、その回のイテレーションでの左辺第１項の微小時間変化を計算するのである。よって、右辺第２項の相互結合項を余裕を持って計算しながら、右辺第１項の利得損失項及び右辺第３項の雑音項を高速に計算することができる。

結合振動子系の計算装置Ｃ３のうち、イテレーション部３１は、利得損失計算部３１１、雑音計算部３１２、相互結合計算部３１３及び振幅位相計算部３１４から構成される。

利得損失計算部３１１は、数２の微分方程式の各イテレーションにおいて、各振動子ｉ毎の利得及び損失に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、利得損失計算部３１１は、数２の微分方程式の各イテレーションにおいて、第ｎ回目のイテレーションの後でのｃ_ｉ ^（ｎ）を入力され、ｃ_ｉ＝Ａｔａｎｈ（ｃ_ｉ／Ｂ）を出力する。ここで、第１回目のイテレーションの前でのｃ_ｉ ^（０）は、０に初期設定される。そして、Ａ及びＢは、数２の右辺第１項のポンプレートｐに応じて設定される。

第１実施形態の利得損失計算部の処理内容を図４に示す。図４の横軸は、利得損失計算部３１１の入力ｃ_ｉｎを示し、図４の縦軸は、利得損失計算部３１１の出力ｃ_ｏｕｔを示す。２つの曲線のうちの一方の曲線は、数２の右辺第１項を計算したときにおける、利得損失計算部３１１の入力ｃ_ｉｎと出力ｃ_ｏｕｔとの間の関係を示す。２つの曲線のうちの他方の曲線は、数２の右辺第１項をｔａｎｈで近似したときにおける、利得損失計算部３１１の入力ｃ_ｉｎと出力ｃ_ｏｕｔとの間の関係を示す。２つの曲線のうち、上記の他方の曲線は、上記の一方の曲線に対して、正当な近似となっていることが分かる。

そこで、利得損失計算部３１１は、図４に示した入力ｃ_ｉｎと出力ｃ_ｏｕｔとの間の関係をルックアップテーブルとして記憶していれば、数２の右辺第１項を効率よく計算することができる。或いは、利得損失計算部３１１は、数２の右辺第１項をルンゲクッタ法又はオイラー丸山法で計算すれば、数２の右辺第１項を精度高く計算することができる。

雑音計算部３１２は、数２の微分方程式の各イテレーションにおいて、各振動子ｉ毎の雑音に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、雑音計算部３１２は、数２の微分方程式の各イテレーションにおいて、利得損失計算部３１１からの出力ｃ_ｉを入力され、ｃ_ｉ＝ｃ_ｉ＋Ｎｏｉｓｅを出力する。ここで、Ｎｏｉｓｅは、数２の右辺第３項に応じて設定される。そして、Ｎｏｉｓｅは、数２の微分方程式の各イテレーションのうち、いずれの段階において付加されてもよい。

相互結合計算部３１３は、数２の微分方程式の各イテレーションにおいて、各振動子ｉの間の相互結合に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、相互結合計算部３１３は、数２の微分方程式の各イテレーションにおいて、雑音計算部３１２からの出力ｃ_ｉを入力され、Σξ_ｉｊｃ_ｊを出力する。ここで、相互結合計算部３１３でのΣξ_ｉｊｃ_ｊの計算時間は、利得損失計算部３１１でのｃ_ｉ＝Ａｔａｎｈ（ｃ_ｉ／Ｂ）の計算時間及び雑音計算部３１２でのｃ_ｉ＝ｃ_ｉ＋Ｎｏｉｓｅの計算時間の合計時間と比べて、イテレーション数回分（例えば、１回分。）だけ長くなっている。

ただし、相互結合計算部３１３は、ξ_ｉｊを表わす行列のうち０となる要素を省略し、ξ_ｉｊを表わす行列を圧縮すれば、数２の右辺第２項を効率よく計算することができる。

振幅位相計算部３１４は、数２の微分方程式の各イテレーションにおいて、利得損失計算部３１１、雑音計算部３１２及び相互結合計算部３１３により各々計算された各振動子ｉ毎の振幅及び位相の微小時間変化に基づいて、各振動子ｉ毎の微小時間後の振幅及び位相を計算する。具体的には、振幅位相計算部３１４は、数２の微分方程式の各イテレーションにおいて、雑音計算部３１２からの出力ｃ_ｉを入力され、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊを入力され、雑音計算部３１２からの出力ｃ_ｉ及び相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊを加算し、第（ｎ＋１）回目のイテレーションの後でのｃ_ｉ ^{（ｎ＋１）}を出力する。ここで、雑音計算部３１２からの出力ｃ_ｉは、第ｎ回目のイテレーションの後でのｃ_ｉ ^（ｎ）を用いて計算されたものであるが、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊは、第（ｎ−ｍ）回目のイテレーションの後でのｃ_ｉ ^{（ｎ−ｍ）}を用いて計算されたものである。ただし、ｍは、イテレーションｍ回分（例えば、１回分。）を表わす。

スピン系計算部３２は、イテレーション部３１により反復されたイテレーションの後に、数２の微分方程式の定常状態における各振動子ｉ毎の位相を計算することにより、結合振動子系の相互結合に対応するイジングモデルの基底状態を計算する。

このように、組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。そして、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。

（第１実施形態：計算結果）
次に、第１実施形態の結合振動子系の計算装置の計算結果を図５に示す。図５の左欄は、完全グラフＫ_４（全振動子数Ｍ＝４）の計算結果を示し、図５の右欄は、ランダムグラフＧ１（全振動子数Ｍ＝８００）の計算結果を示す。図５の横軸は、イテレーション回数を示し、図５の縦軸は、振幅位相計算部３１４の出力ｃ_ｉ ^（ｎ）を示す。

図５の左欄の完全グラフＫ_４では、第１回目のイテレーションの前でのｃ_ｉ ^（０）を０に初期設定し、数２の右辺第１項のポンプレートをｐ＝２．０に設定し、利得損失計算部３１１からの出力をｃ_ｏｕｔ＝２．２ｔａｎｈ（ｃ_ｉｎ／１．８）に設定した。すると、イテレーション回数〜３０回目において、ｃ_ｉ ^（ｎ）が０から立ち上がり、イテレーション回数〜５０回目において、ｃ_ｉ ^（ｎ）が定常状態に落ち着き、完全グラフＫ_４の最大のカット数が得られ、完全グラフＫ_４に対応するイジングモデルの基底状態が得られた。

図５の右欄のランダムグラフＧ１では、第１回目のイテレーションの前でのｃ_ｉ ^（０）を０に初期設定し、数２の右辺第１項のポンプレートをｐ＝１．１に設定し、利得損失計算部３１１からの出力をｃ_ｏｕｔ＝２．３５ｔａｎｈ（ｃ_ｉｎ／２．４６）に設定した。すると、イテレーション回数〜１０回目において、ｃ_ｉ ^（ｎ）が０から立ち上がり、イテレーション回数〜１００回目において、ｃ_ｉ ^（ｎ）が定常状態に落ち着き、ランダムグラフＧ１の最大のカット数が得られ、ランダムグラフＧ１に対応するイジングモデルの基底状態が得られた。

ここで、雑音計算部３１２は、数２の微分方程式の初期状態のみにおいて、各振動子ｉ毎の雑音に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算してもよい。具体的には、雑音計算部３１２は、図５の左欄の完全グラフＫ_４では、イテレーション回数〜３０回目において、計算を停止してもよく、図５の右欄のランダムグラフＧ１では、イテレーション回数〜１０回目において、計算を停止してもよい。

よって、数２の微分方程式の定常状態として、複数の結合振動子の振幅が０である状態を回避することができる。そして、複数の結合振動子の振幅が０である初期状態から脱出することができた後には、雑音計算部３１２の計算負担をなくすことができる。

さらに、イテレーション部３１は、前回のイテレーションにおける各振動子ｉ毎の微小時間後の振幅及び位相が、数２の微分方程式の定常状態に所定の精度で接近したときに、次回のイテレーションにおける各振動子ｉ毎の微小時間後の振幅及び位相の計算を停止してもよい。具体的には、イテレーション部３１は、図５の左欄の完全グラフＫ_４では、イテレーション回数〜５０回目において、計算を停止してもよく、図５の右欄のランダムグラフＧ１では、イテレーション回数〜１００回目において、計算を停止してもよい。

より一般的には、イテレーション部３１は、所定の割合（例えば、全振動子の９０％。）の振動子について、ｃ_ｉ ^（ｎ）が所定の閾値（例えば、飽和振幅の９０％。）を超えたときに、計算を停止してもよい。よって、複数の結合振動子の振幅が有限値である定常状態に接近することができた後には、イテレーション部３１の計算負担をなくすことができる。

（第１実施形態：並列計算）
次に、第１実施形態のイテレーション部の第１構成を図６に示す。イテレーション部３１は、複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍ及びＦＩＦＯ３３から構成される。複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍ及びＦＩＦＯ３３は、１プロセッサ又は１ＦＰＧＡに収容される。

複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍは、数２の微分方程式の各イテレーションにおいて、各振動子１、・・・、ｍ、・・・、Ｍ毎の並列処理により、各振動子１、・・・、ｍ、・・・、Ｍ毎の微小時間後の振幅及び位相を計算する。

図６の以下の説明では、イテレーション部３１−ｍについて説明する。イテレーション部３１−ｍは、利得損失計算部３１１−ｍ、雑音計算部３１２−ｍ、相互結合計算部３１３−ｍ及び振幅位相計算部３１４−ｍから構成される。

利得損失計算部３１１−ｍは、ＦＩＦＯ３３（詳細については後述する。）に格納されたｃ_１ ^（ｎ）、・・・、ｃ_ｍ ^（ｎ）、・・・、ｃ_Ｍ ^（ｎ）をこの順序で入力され、ｃ_１＝Ａｔａｎｈ（ｃ_１／Ｂ）、・・・、ｃ_ｍ＝Ａｔａｎｈ（ｃ_ｍ／Ｂ）、・・・、ｃ_Ｍ＝Ａｔａｎｈ（ｃ_Ｍ／Ｂ）をこの順序で出力する。ここで、利得損失計算部３１１−ｍは、複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍの間で共有されていてもよい。

雑音計算部３１２−ｍは、利得損失計算部３１１−ｍからの出力ｃ_１、・・・、ｃ_ｍ、・・・、ｃ_Ｍをこの順序で入力され、ｃ_１＝ｃ_１＋Ｎｏｉｓｅ、・・・、ｃ_ｍ＝ｃ_ｍ＋Ｎｏｉｓｅ、・・・、ｃ_Ｍ＝ｃ_Ｍ＋Ｎｏｉｓｅをこの順序で出力する。ここで、雑音計算部３１２−ｍは、複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍの間で共有されていてもよい。

相互結合計算部３１３−ｍは、雑音計算部３１２−ｍからの出力ｃ_１、・・・、ｃ_ｍ、・・・、ｃ_Ｍをこの順序で入力され、Σξ_ｍｊｃ_ｊを出力する。ここで、相互結合計算部３１３−ｍは、ＦＩＦＯ３１５−ｍ、乗算器３１６−ｍ及び加算器３１７−ｍから構成される。

乗算器３１６−ｍは、雑音計算部３１２−ｍからの出力ｃ_１、・・・、ｃ_ｍ、・・・、ｃ_Ｍをこの順序で入力され、ＦＩＦＯ３１５−ｍに格納されたξ_ｍ１、・・・、ξ_ｍｍ（＝０）、・・・、ξ_ｍＭをこの順序で入力され、ξ_ｍ１ｃ_１、・・・、ξ_ｍｍｃ_ｍ（＝０）、・・・、ξ_ｍＭｃ_Ｍをこの順序で出力する。ここで、乗算器３１６−ｍでのξ_ｍ１ｃ_１、・・・、ξ_ｍｍｃ_ｍ（＝０）、・・・、ξ_ｍＭｃ_Ｍの計算時間は、利得損失計算部３１１−ｍでのｃ_ｍ＝Ａｔａｎｈ（ｃ_ｍ／Ｂ）の計算時間及び雑音計算部３１２−ｍでのｃ_ｍ＝ｃ_ｍ＋Ｎｏｉｓｅの計算時間の合計時間と比べて、イテレーション数回分（例えば、１回分。）だけ長くなっている。

加算器３１７−ｍは、乗算器３１６−ｍからの出力ξ_ｍ１ｃ_１、・・・、ξ_ｍｍｃ_ｍ（＝０）、・・・、ξ_ｍＭｃ_Ｍをこの順序で入力され、Σξ_ｍｊｃ_ｊを出力する。

振幅位相計算部３１４−ｍは、雑音計算部３１２−ｍからの出力ｃ_ｍを入力され、加算器３１７−ｍからの出力Σξ_ｍｊｃ_ｊを入力され、雑音計算部３１２−ｍからの出力ｃ_ｍ及び加算器３１７−ｍからの出力Σξ_ｍｊｃ_ｊを加算し、ｃ_ｍ ^{（ｎ＋１）}を出力する。ここで、雑音計算部３１２−ｍからの出力ｃ_ｍは、第ｎ回目のイテレーションの後でのｃ_ｍ ^（ｎ）を用いて計算されたものであるが、加算器３１７−ｍからの出力Σξ_ｍｊｃ_ｊは、第（ｎ−ｍ）回目のイテレーションの後でのｃ_ｍ ^{（ｎ−ｍ）}を用いて計算されたものである。ただし、ｍは、イテレーションｍ回分（例えば、１回分。）を表わす。

最初に、イテレーション部３１−１は、ｃ_１ ^{（ｎ＋１）}をＦＩＦＯ３３に格納する。・・・次に、イテレーション部３１−ｍは、ｃ_ｍ ^{（ｎ＋１）}をＦＩＦＯ３３に格納する。・・・最後に、イテレーション部３１−Ｍは、ｃ_Ｍ ^{（ｎ＋１）}をＦＩＦＯ３３に格納する。

ここで、結合振動子系の全振動子数をＭとすると、各振動子毎の並列処理を行わないときには、２体相互結合の計算量は、Ｍ^２のオーダーに増える。しかし、各振動子毎の並列処理を行なうときには、２体相互結合の計算量は、ｌｏｇ_２Ｍのオーダーに減らすことができる。

第１実施形態のイテレーション部の第２構成を図７に示す。イテレーション部３１は、複数のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩから構成される。各々のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩは、１プロセッサ又は１ＦＰＧＡに収容される。

図６に示した第１構成に加えて、各々のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩは、数２の微分方程式の各イテレーションにおいて、結合振動子系のうちの各々に特定のグループに属する振動子の微小時間後の振幅及び位相のデータを、その他のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩとの間で共有する。

イテレーション部３１−Ｉは、ＦＩＦＯ３３−Ｉ、イテレーション部３４−Ｉ及びＦＩＦＯ３５−Ｉから構成される。イテレーション部３４−Ｉは、ｃ_１、・・・、ｃ_３Ｍのうち、ｃ_１、ｃ_２、・・・、ｃ_Ｍについて、図６に示した複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍと同様の処理を実行する。

イテレーション部３４−Ｉは、ＦＩＦＯ３３−Ｉ（詳細については後述する。）に格納されたｃ_１ ^（ｎ）、・・・、ｃ_Ｍ ^（ｎ）、ｃ_Ｍ＋１ ^（ｎ）、・・・、ｃ_２Ｍ ^（ｎ）、ｃ_２Ｍ＋１ ^（ｎ）、・・・、ｃ_３Ｍ ^（ｎ）をこの順序で入力され、ｃ_１ ^{（ｎ＋１）}、ｃ_２ ^{（ｎ＋１）}、・・・、ｃ_Ｍ ^{（ｎ＋１）}をこの順序で出力し、ＦＩＦＯ３３−Ｉに格納する前にＦＩＦＯ３５−Ｉに格納する。

イテレーション部３１−ＩＩは、ＦＩＦＯ３３−ＩＩ、イテレーション部３４−ＩＩ及びＦＩＦＯ３５−ＩＩから構成される。イテレーション部３４−ＩＩは、ｃ_１、・・・、ｃ_３Ｍのうち、ｃ_Ｍ＋１、ｃ_Ｍ＋２、・・・、ｃ_２Ｍについて、図６に示した複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍと同様の処理を実行する。

イテレーション部３４−ＩＩは、ＦＩＦＯ３３−ＩＩ（詳細については後述する。）に格納されたｃ_１ ^（ｎ）、・・・、ｃ_Ｍ ^（ｎ）、ｃ_Ｍ＋１ ^（ｎ）、・・・、ｃ_２Ｍ ^（ｎ）、ｃ_２Ｍ＋１ ^（ｎ）、・・・、ｃ_３Ｍ ^（ｎ）をこの順序で入力され、ｃ_Ｍ＋１ ^{（ｎ＋１）}、ｃ_Ｍ＋２ ^{（ｎ＋１）}、・・・、ｃ_２Ｍ ^{（ｎ＋１）}をこの順序で出力し、ＦＩＦＯ３３−ＩＩに格納する前にＦＩＦＯ３５−ＩＩに格納する。

イテレーション部３１−ＩＩＩは、ＦＩＦＯ３３−ＩＩＩ、イテレーション部３４−ＩＩＩ及びＦＩＦＯ３５−ＩＩＩから構成される。イテレーション部３４−ＩＩＩは、ｃ_１、・・・、ｃ_３Ｍのうち、ｃ_２Ｍ＋１、ｃ_２Ｍ＋２、・・・、ｃ_３Ｍについて、図６に示した複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍと同様の処理を実行する。

イテレーション部３４−ＩＩＩは、ＦＩＦＯ３３−ＩＩＩ（詳細については後述する。）に格納されたｃ_１ ^（ｎ）、・・・、ｃ_Ｍ ^（ｎ）、ｃ_Ｍ＋１ ^（ｎ）、・・・、ｃ_２Ｍ ^（ｎ）、ｃ_２Ｍ＋１ ^（ｎ）、・・・、ｃ_３Ｍ ^（ｎ）をこの順序で入力され、ｃ_２Ｍ＋１ ^{（ｎ＋１）}、ｃ_２Ｍ＋２ ^{（ｎ＋１）}、・・・、ｃ_３Ｍ ^{（ｎ＋１）}をこの順序で出力し、ＦＩＦＯ３３−ＩＩＩに格納する前にＦＩＦＯ３５−ＩＩＩに格納する。

最初に、ｃ_１ ^{（ｎ＋１）}、ｃ_２ ^{（ｎ＋１）}、・・・、ｃ_Ｍ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。次に、ｃ_Ｍ＋１ ^{（ｎ＋１）}、ｃ_Ｍ＋２ ^{（ｎ＋１）}、・・・、ｃ_２Ｍ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。最後に、ｃ_２Ｍ＋１ ^{（ｎ＋１）}、ｃ_２Ｍ＋２ ^{（ｎ＋１）}、・・・、ｃ_３Ｍ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩＩＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。ここで、ＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩは、単一のＦＩＦＯであってもよい。

このように、結合振動子系の全振動子数が膨大な数になるときでも、各々のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩを各々のプロセッサに実装することができる。

第１実施形態のイテレーション部の第３構成を図８に示す。イテレーション部３１は、複数のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩから構成される。各々のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩは、１プロセッサ又は１ＦＰＧＡに収容される。

図７に示した第２構成に加えて、前回のイテレーションにおける上記のデータ（ｃ_１ ^{（ｎ＋１）}、・・・ｃ_３Ｍ ^{（ｎ＋１）}）の共有処理が、次回のイテレーションにおける計算処理を停止／中断させないように、各々のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩは、結合振動子系のうちの各々に特定のグループに属する振動子を割り当てられる。

イテレーション部３１−Ｉは、ＦＩＦＯ３３−Ｉ、イテレーション部３４−Ｉ及びＦＩＦＯ３５−Ｉから構成される。イテレーション部３４−Ｉは、ｃ_１、・・・、ｃ_３Ｍのうち、ｃ_１、ｃ_４、・・・、ｃ_３Ｍ−２について、図６に示した複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍと同様の処理を実行する。

イテレーション部３４−Ｉは、ＦＩＦＯ３３−Ｉ（詳細については後述する。）に格納されたｃ_１ ^（ｎ）、・・・、ｃ_Ｍ ^（ｎ）、ｃ_Ｍ＋１ ^（ｎ）、・・・、ｃ_２Ｍ ^（ｎ）、ｃ_２Ｍ＋１ ^（ｎ）、・・・、ｃ_３Ｍ ^（ｎ）をこの順序で入力され、ｃ_１ ^{（ｎ＋１）}、ｃ_４ ^{（ｎ＋１）}、・・・、ｃ_３Ｍ−２ ^{（ｎ＋１）}をこの順序で出力し、ＦＩＦＯ３３−Ｉに格納する前にＦＩＦＯ３５−Ｉに格納する。

イテレーション部３１−ＩＩは、ＦＩＦＯ３３−ＩＩ、イテレーション部３４−ＩＩ及びＦＩＦＯ３５−ＩＩから構成される。イテレーション部３４−ＩＩは、ｃ_１、・・・、ｃ_３Ｍのうち、ｃ_２、ｃ_５、・・・、ｃ_３Ｍ−１について、図６に示した複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍと同様の処理を実行する。

イテレーション部３４−ＩＩは、ＦＩＦＯ３３−ＩＩ（詳細については後述する。）に格納されたｃ_１ ^（ｎ）、・・・、ｃ_Ｍ ^（ｎ）、ｃ_Ｍ＋１ ^（ｎ）、・・・、ｃ_２Ｍ ^（ｎ）、ｃ_２Ｍ＋１ ^（ｎ）、・・・、ｃ_３Ｍ ^（ｎ）をこの順序で入力され、ｃ_２ ^{（ｎ＋１）}、ｃ_５ ^{（ｎ＋１）}、・・・、ｃ_３Ｍ−１ ^{（ｎ＋１）}をこの順序で出力し、ＦＩＦＯ３３−ＩＩに格納する前にＦＩＦＯ３５−ＩＩに格納する。

イテレーション部３１−ＩＩＩは、ＦＩＦＯ３３−ＩＩＩ、イテレーション部３４−ＩＩＩ及びＦＩＦＯ３５−ＩＩＩから構成される。イテレーション部３４−ＩＩＩは、ｃ_１、・・・、ｃ_３Ｍのうち、ｃ_３、ｃ_６、・・・、ｃ_３Ｍについて、図６に示した複数のイテレーション部３１−１、・・・、３１−ｍ、・・・、３１−Ｍと同様の処理を実行する。

イテレーション部３４−ＩＩＩは、ＦＩＦＯ３３−ＩＩＩ（詳細については後述する。）に格納されたｃ_１ ^（ｎ）、・・・、ｃ_Ｍ ^（ｎ）、ｃ_Ｍ＋１ ^（ｎ）、・・・、ｃ_２Ｍ ^（ｎ）、ｃ_２Ｍ＋１ ^（ｎ）、・・・、ｃ_３Ｍ ^（ｎ）をこの順序で入力され、ｃ_３ ^{（ｎ＋１）}、ｃ_６ ^{（ｎ＋１）}、・・・、ｃ_３Ｍ ^{（ｎ＋１）}をこの順序で出力し、ＦＩＦＯ３３−ＩＩＩに格納する前にＦＩＦＯ３５−ＩＩＩに格納する。

最初に、ｃ_１ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。次に、ｃ_２ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。次に、ｃ_３ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩＩＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。・・・次に、ｃ_３Ｍ−２ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。次に、ｃ_３Ｍ−１ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。最後に、ｃ_３Ｍ ^{（ｎ＋１）}が、ＦＩＦＯ３５−ＩＩＩからＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩへと転送される。ここで、ＦＩＦＯ３３−Ｉ、３３−ＩＩ、３３−ＩＩＩは、単一のＦＩＦＯであってもよい。

このように、各々のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩを各々のプロセッサに実装するときでも、各々のイテレーション部３１−Ｉ、３１−ＩＩ、３１−ＩＩＩの間の上記のデータ（ｃ_１ ^{（ｎ＋１）}、・・・ｃ_３Ｍ ^{（ｎ＋１）}）の共有処理を円滑に実行することができる。

（第１実施形態：ビット量子化）
次に、第１実施形態の振幅位相計算部の量子化を図９に示す。図９の左欄は、振幅位相計算部３１４の量子化の第１形態を示し、図９の右欄は、振幅位相計算部３１４の量子化の第２形態を示す。図９の横軸は、振幅位相計算部３１４の量子化前の入力を示し、図９の縦軸は、振幅位相計算部３１４の量子化後の出力を示す。

振幅位相計算部３１４は、数２の微分方程式の各イテレーションにおいて、所定の精度で量子化したうえで、各振動子ｉ毎の微小時間後の振幅及び位相を計算する。具体的には、振幅位相計算部３１４は、数２の微分方程式の各イテレーションにおいて、雑音計算部３１２からの出力ｃ_ｉ及び相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊの加算結果を所定の精度で量子化し、第（ｎ＋１）回目のイテレーションの後でのｃ_ｉ ^{（ｎ＋１）}を出力する。

図９の左欄に示した第１形態では、振幅位相計算部３１４の出力ビット数ｘは、４である。そして、振幅位相計算部３１４の出力は、０を含まず、振幅位相計算部３１４の出力の階調数は、２^ｘ＝２^４＝１６である。よって、数２の微分方程式の定常状態として、複数の結合振動子の振幅が０である状態を回避することができる。

図９の右欄に示した第２形態でも、振幅位相計算部３１４の出力ビット数ｘは、４である。しかし、振幅位相計算部３１４の出力は、０を含み、振幅位相計算部３１４の出力の階調数は、２^ｘ−１＝２^４−１＝１５である。よって、数２の微分方程式の初期状態において、結合振動子の位相が正答の位相からずれた位相になる状態を回避することができる（図５の右欄に示したｃ_ｉ ^（ｎ）の正／負から負／正への変化を参照。）。

このように、結合振動子系の計算装置Ｃ３の計算負担を減らすことができるとともに、結合振動子系の計算装置Ｃ３のメモリ（ＦＩＦＯ）容量を減らすことができる。

振幅位相計算部３１４により計算される各振動子ｉ毎の微小時間後の振幅についての量子化カットオフは、振幅位相計算部３１４により計算される各振動子ｉ毎の微小時間後の振幅についての量子化の上記の所定の精度を満足するように設定されている。具体的には、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化カットオフは、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化の上記の所定の精度を満足する。

図９の左欄及び右欄に示した第１、２形態では、振幅位相計算部３１４の量子化カットオフ±Ｎは、±２である。よって、振幅位相計算部３１４の量子化精度は、Ｎ×２／２^ｘ＝２×２／２^４＝０．２５程度である。もっとも、図１０〜１２に示すように、振幅位相計算部３１４の量子化精度は、０．２５程度では正答を得るには低過ぎる。

振幅位相計算部３１４により計算される各振動子ｉ毎の微小時間後の振幅についての量子化カットオフは、数２の微分方程式の定常状態における各振動子ｉ毎の振幅と比べて大きく設定されている。具体的には、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化カットオフは、振幅位相計算部３１４からの定常出力ｃ_ｉ ^{（ｎ＋１）}と比べて大きい。

図９の左欄及び右欄に示した第１、２形態では、振幅位相計算部３１４の量子化カットオフ±Ｎは、±２である。そして、図５に示したように、振幅位相計算部３１４からの定常出力ｃ_ｉ ^{（ｎ＋１）}は、±１〜２程度である。よって、振幅位相計算部３１４の量子化カットオフＮは、振幅位相計算部３１４からの定常出力ｃ_ｉ ^{（ｎ＋１）}と比べて大きい。

このように、各振動子毎の振幅の分解能を低くし過ぎず、かつ、各振動子毎の振幅の制限を厳しくし過ぎず、結合振動子系の計算装置Ｃ３の計算負担を減らすことができるとともに、結合振動子系の計算装置Ｃ３のメモリ（ＦＩＦＯ）容量を減らすことができる。

相互結合計算部３１３により計算される各振動子ｉ毎の振幅の微小時間変化についての量子化ビット数は、振幅位相計算部３１４により計算される各振動子ｉ毎の微小時間後の振幅についての量子化ビット数と比べて、少なくともｌｏｇ_２Ｍだけ大きく設定されている。ただし、Ｍは、結合振動子系の全振動子数である。具体的には、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊについての量子化ビット数は、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化ビット数と比べて、少なくともｌｏｇ_２Ｍだけ大きい。

仮に、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊについての量子化ビット数を、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化ビット数と比べて、同等にする。すると、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊは、Ｍ個のｃ_ｊを重み付け加算したものであるため、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊについての量子化精度は、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化精度と比べて、低くなる。

そこで、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊについての量子化ビット数を、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化ビット数と比べて、少なくともｌｏｇ_２Ｍだけ大きくする。すると、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊは、Ｍ個のｃ_ｊを重み付け加算したものであるところ、相互結合計算部３１３からの出力Σξ_ｉｊｃ_ｊについての量子化精度は、振幅位相計算部３１４からの出力ｃ_ｉ ^{（ｎ＋１）}についての量子化精度と比べて、少なくとも同等になる。

このように、各振動子毎の振幅の分解能と比べて、各振動子毎の相互結合の大きさの分解能を低くし過ぎず、結合振動子系の計算装置Ｃ３の計算負担を減らすことができるとともに、結合振動子系の計算装置Ｃ３のメモリ（ＦＩＦＯ）容量を減らすことができる。

第１実施形態の結合振動子系の計算装置における、完全グラフＫ_４（全振動子数Ｍ＝４）の計算結果を図１０に示す。１００００回の試行のうち、正答率が高かったときほど、塗りつぶしを薄くし、正答率が低かったときほど、塗りつぶしを濃くする。

図１０の左欄の横軸は、振幅位相計算部３１４の出力ビット数ｘを示し、図１０の左欄の縦軸は、利得損失計算部３１１の出力ビット数ｙを示す。ただし、相互結合計算部３１３の出力ビット数ｚは、３２に固定されており、振幅位相計算部３１４の量子化カットオフＮは、２．０に固定されている。ｘ及びｙがともに小さいときや、ｘが大きくてもｙが小さいときや、ｙが大きくてもｘが小さいときは、正答率が低くなった。ｘ及びｙがほぼ等しくともに大きい（例えば、１５程度。）ときは、正答率が高くなった。

図１０の中欄の横軸は、振幅位相計算部３１４の出力ビット数ｘ＝利得損失計算部３１１の出力ビット数ｙを示し、図１０の中欄の縦軸は、相互結合計算部３１３の出力ビット数ｚを示す。ただし、振幅位相計算部３１４の量子化カットオフＮは、２．０に固定されている。ｘ＝ｙが大きくなるにつれて、ｚがより大きくならなければ、相互結合が正確に計算されず正答率が低くなった。ｘ＝ｙが大きい（例えば、１５程度。）うえで、ｚがｘ＝ｙにほぼ等しいときは、相互結合がやや正確に計算され正答率がやや高くなった。ｘ＝ｙが大きい（例えば、１５程度。）うえで、ｚがｘ＝ｙよりｌｏｇ_２Ｍ＝ｌｏｇ_２４＝２だけ大きいときは、相互結合がより正確に計算され正答率がより高くなった。

図１０の右欄の横軸は、振幅位相計算部３１４の出力ビット数ｘを示し、図１０の右欄の縦軸は、振幅位相計算部３１４の量子化カットオフｌｏｇ_１０Ｎを示す。ただし、利得損失計算部３１１の出力ビット数ｙは、振幅位相計算部３１４の出力ビット数ｘに等しく、相互結合計算部３１３の出力ビット数ｚは、３２に固定されている。ｘが小さくなるにつれて、Ｎがより小さくならなければ、量子化精度が低くなり正答率が低くなった。ｘが図１０の右欄に示したいずれの値であっても、Ｎが小さすぎるときは、振動子振幅が制限され正答率が低くなった。ｘが大きい（例えば、１５程度。）うえで、Ｎが１程度であるときは、量子化精度が高くなり振動子振幅が制限されず、正答率がより高くなった。

第１実施形態の結合振動子系の計算装置における、メビウスラダーＭ_１６（全振動子数Ｍ＝１６）の計算結果を図１１に示す。１００００回の試行のうち、正答率が高かったときほど、塗りつぶしを薄くし、正答率が低かったときほど、塗りつぶしを濃くする。

図１１の左欄の横軸は、振幅位相計算部３１４の出力ビット数ｘを示し、図１１の左欄の縦軸は、利得損失計算部３１１の出力ビット数ｙを示す。ただし、相互結合計算部３１３の出力ビット数ｚは、３２に固定されており、振幅位相計算部３１４の量子化カットオフＮは、２．０に固定されている。ｘ及びｙがともに小さいときや、ｘが大きくてもｙが小さいときや、ｙが大きくてもｘが小さいときは、正答率が低くなった。ｘ及びｙがほぼ等しくともに大きい（例えば、１５程度。）ときは、正答率が高くなった。

図１１の中欄の横軸は、振幅位相計算部３１４の出力ビット数ｘ＝利得損失計算部３１１の出力ビット数ｙを示し、図１１の中欄の縦軸は、相互結合計算部３１３の出力ビット数ｚを示す。ただし、振幅位相計算部３１４の量子化カットオフＮは、２．０に固定されている。ｘ＝ｙが大きくなるにつれて、ｚがより大きくならなければ、相互結合が正確に計算されず正答率が低くなった。ｘ＝ｙが大きい（例えば、１５程度。）うえで、ｚがｘ＝ｙにほぼ等しいときは、相互結合がやや正確に計算され正答率がやや高くなった。ｘ＝ｙが大きい（例えば、１５程度。）うえで、ｚがｘ＝ｙよりｌｏｇ_２Ｍ＝ｌｏｇ_２１６＝４だけ大きいときは、相互結合がより正確に計算され正答率がより高くなった。

図１１の右欄の横軸は、振幅位相計算部３１４の出力ビット数ｘを示し、図１１の右欄の縦軸は、振幅位相計算部３１４の量子化カットオフｌｏｇ_１０Ｎを示す。ただし、利得損失計算部３１１の出力ビット数ｙは、振幅位相計算部３１４の出力ビット数ｘに等しく、相互結合計算部３１３の出力ビット数ｚは、３２に固定されている。ｘが小さくなるにつれて、Ｎがより小さくならなければ、量子化精度が低くなり正答率が低くなった。ｘが図１１の右欄に示したいずれの値であっても、Ｎが小さすぎるときは、振動子振幅が制限され正答率が低くなった。ｘが大きい（例えば、１５程度。）うえで、Ｎが１程度であるときは、量子化精度が高くなり振動子振幅が制限されず、正答率がより高くなった。

第１実施形態の結合振動子系の計算装置における、ランダムグラフＧ１（全振動子数Ｍ＝８００）の計算結果を図１２に示す。１０００回の試行のうち、カット数が高かったときほど、塗りつぶしを薄くし、カット数が低かったときほど、塗りつぶしを濃くする。

図１２の左欄の横軸は、振幅位相計算部３１４の出力ビット数ｘを示し、図１２の左欄の縦軸は、利得損失計算部３１１の出力ビット数ｙを示す。ただし、相互結合計算部３１３の出力ビット数ｚは、１６に固定されており、振幅位相計算部３１４の量子化カットオフＮは、１．０に固定されている。ｘ及びｙがともに小さいときや、ｘが大きくてもｙが小さいときや、ｙが大きくてもｘが小さいときは、カット数が低くなった。ｘ及びｙがほぼ等しくともに大きい（例えば、１５程度。）ときは、カット数が高くなった。

図１２の中欄の横軸は、振幅位相計算部３１４の出力ビット数ｘ＝利得損失計算部３１１の出力ビット数ｙを示し、図１２の中欄の縦軸は、相互結合計算部３１３の出力ビット数ｚを示す。ただし、振幅位相計算部３１４の量子化カットオフＮは、１．０に固定されている。ｘ＝ｙが大きくなるにつれて、ｚがより大きくならなければ、相互結合が正確に計算されずカット数が低くなった。ｘ＝ｙが大きい（例えば、１５程度。）うえで、ｚがｘ＝ｙにほぼ等しいときは、相互結合がやや正確に計算されカット数がやや高くなった。ｘ＝ｙが大きい（例えば、１５程度。）うえで、ｚがｘ＝ｙよりｌｏｇ_２Ｍ＝ｌｏｇ_２８００〜１０だけ大きいときは、相互結合がより正確に計算されカット数がより高くなった。

図１２の右欄の横軸は、振幅位相計算部３１４の出力ビット数ｘを示し、図１２の右欄の縦軸は、振幅位相計算部３１４の量子化カットオフｌｏｇ_１０Ｎを示す。ただし、利得損失計算部３１１の出力ビット数ｙは、振幅位相計算部３１４の出力ビット数ｘに等しく、相互結合計算部３１３の出力ビット数ｚは、３２に固定されている。ｘが小さくなるにつれて、Ｎがより小さくならなければ、量子化精度が低くなりカット数が低くなった。ｘが図１２の右欄に示したいずれの値であっても、Ｎが小さすぎるときは、振動子振幅が制限されカット数が低くなった。ｘが大きい（例えば、１５程度。）うえで、Ｎが１程度であるときは、量子化精度が高くなり振動子振幅が制限されず、カット数がより高くなった。

（第２実施形態：計算方法）
第２実施形態の結合振動子系の計算装置の構成を図１３に示す。結合振動子系の計算装置Ｃ４は、イテレーション部４１及びスピン系計算部４２から構成される。結合振動子系の計算装置Ｃ４に対して、結合振動子系の計算プログラムをインストールすることにより、イテレーション部４１及びスピン系計算部４２を実装することができる。

イテレーション部４１は、結合振動子系のダイナミクスを記述する微分方程式（数２）の各イテレーションにおいて、各振動子ｉ毎の微小時間後の振幅及び位相（数２のｃ_ｉの絶対値及び符号）を計算する。よって、イテレーション部４１は、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算することができる。

スピン系計算部４２は、イテレーション部４１により反復されたイテレーションの後に、微分方程式（数２）の定常状態における各振動子ｉ毎の位相（数２のｃ_ｉの符号）を計算することにより、結合振動子系の相互結合（数２のξ_ｉｊ）に対応するイジングモデル（数１）の基底状態を計算する。よって、スピン系計算部４２は、イジングモデルの基底状態を計算することにより、組み合わせ最適化問題の最適解を求めることができる。

ここで、結合振動子系の全振動子数をＭとすると、数２の微分方程式のうち、右辺第１項の利得損失項及び右辺第３項の雑音項の計算量は、Ｍに比例するが、右辺第２項の相互結合項の計算量は、Ｍ^２に比例する。しかし、数２の微分方程式のうち、右辺第２項の相互結合項の計算時間が、右辺第１項の利得損失項及び右辺第３項の雑音項の計算時間と比べて、同程度であるときには、図１３に示したように各イテレーションを行なえばよい。

つまり、ある回のイテレーションの後でのｃ_ｉを用いて計算した、右辺第１項の利得損失項、右辺第２項の相互結合項及び右辺第３項の雑音項を加算することにより、その回のイテレーションでの左辺第１項の微小時間変化を計算するのである。よって、右辺第２項の相互結合項を図３のような遅延（例えば、イテレーション数回分。）を伴わず計算しながら、右辺第１項の利得損失項及び右辺第３項の雑音項を高速に計算することができる。

結合振動子系の計算装置Ｃ４のうち、イテレーション部４１は、利得損失計算部４１１、相互結合計算部４１２、雑音計算部４１３及び振幅位相計算部４１４から構成される。

利得損失計算部４１１は、数２の微分方程式の各イテレーションにおいて、各振動子ｉ毎の利得及び損失に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、利得損失計算部４１１は、数２の微分方程式の各イテレーションにおいて、第ｎ回目のイテレーションの後でのｃ_ｉ ^（ｎ）を入力され、数２の右辺第１項を出力する。ここで、第１回目のイテレーションの前でのｃ_ｉ ^（０）は、０に初期設定される。

雑音計算部４１３は、数２の微分方程式の各イテレーションにおいて、各振動子ｉ毎の雑音に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、雑音計算部４１３は、数２の微分方程式の各イテレーションにおいて、第ｎ回目のイテレーションの後でのｃ_ｉ ^（ｎ）を入力され、数２の右辺第３項を出力する。

相互結合計算部４１２は、数２の微分方程式の各イテレーションにおいて、各振動子ｉの間の相互結合に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、相互結合計算部４１２は、数２の微分方程式の各イテレーションにおいて、第ｎ回目のイテレーションの後でのｃ_ｉ ^（ｎ）を入力され、数２の右辺第２項を出力する。ここで、相互結合計算部４１２での数２の右辺第２項の計算時間は、利得損失計算部４１１及び雑音計算部４１３での数２の右辺第１、３項の計算時間と比べて、同程度である。

振幅位相計算部４１４は、数２の微分方程式の各イテレーションにおいて、利得損失計算部４１１、雑音計算部４１３及び相互結合計算部４１２により各々計算された各振動子ｉ毎の振幅及び位相の微小時間変化に基づいて、各振動子ｉ毎の微小時間後の振幅及び位相を計算する。具体的には、振幅位相計算部４１４は、数２の微分方程式の各イテレーションにおいて、利得損失計算部４１１、雑音計算部４１３及び相互結合計算部４１２からの数２の右辺第１、３、２項を入力され、利得損失計算部４１１、雑音計算部４１３及び相互結合計算部４１２からの数２の右辺第１、３、２項を加算し、第（ｎ＋１）回目のイテレーションの後でのｃ_ｉ ^{（ｎ＋１）}を出力する。ここで、利得損失計算部４１１、雑音計算部４１３及び相互結合計算部４１２からの数２の右辺第１、３、２項は、第ｎ回目のイテレーションの後でのｃ_ｉ ^（ｎ）を用いて計算されたものである。

スピン系計算部４２は、イテレーション部４１により反復されたイテレーションの後に、数２の微分方程式の定常状態における各振動子ｉ毎の位相を計算することにより、結合振動子系の相互結合に対応するイジングモデルの基底状態を計算する。

このように、組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。そして、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。

第２実施形態においても、第１実施形態の図６〜８と同様に、並列計算を行なうことができ、第１実施形態の図９〜１２と同様に、ビット量子化を行なうことができる。

（第３実施形態：計算原理）
第３実施形態では、機械学習のボルツマンサンプリングを行なうために、機械学習をＸＹモデルへマッピングし、ＸＹモデルの平衡状態を計算する。または、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算する。

そこで、複数の振動子で複数のＸＹ面内スピンを実装する。ここで、複数の振動子の間の相互結合を行なうことにより、複数のＸＹ面内スピンの間の相互結合を実装する。そして、複数の振動子全体の様々な発振モードのサンプリングを行なうことにより、複数のＸＹ面内スピン全体の様々なエネルギーのサンプリングを行なう。

まず、複数の振動子全体の様々な発振モードのサンプリングを行なうことにより、複数のＸＹ面内スピン全体の様々なエネルギーのサンプリングを行なうことができる理由について説明する。ＸＹモデルのハミルトニアンは、数７で表される。ここで、Ｊ_ｉｊは、スピンｉ、ｊの間の相互結合係数である。そして、θ_ｉ、θ_ｊは、それぞれ、スピンｉ、ｊのＸＹ面内での連続値をとる位相（０ｒａｄから２πｒａｄまで）を表わす。

複数の振動子の時間発展を表わす複素ランジュバン方程式は、数８により表わされる。ここで、ｄＡ_ｉは、振動子ｉの信号成分についての複素振幅変化を表わし、ｄｔは、時間変化を表わし、ｄＷ_ｉは、振動子ｉの雑音成分についての複素ウィーナー過程を表わす。

ｇ（Ａ_ｉ）は、振動子ｉの強度増幅率であり、数９により表わされる。ここで、ｇ_０は、小信号利得を表わし、ｎ_０は、飽和強度を表わす。よって、（１／２）ｇ（Ａ_ｉ）Ａ_ｉｄｔは、振動子ｉの信号成分について、利得に起因する複素振幅変化を表わす。

γ_ｃは、振動子の強度減衰率であり、数１０により表わされる。ここで、ωは、振動子の発振周波数を表わし、Ｑは、振動子の共振Ｑ値を表わす。ただし、ω及びＱは、すべての振動子について、同一の値をとる。よって、−（１／２）γ_ｃＡ_ｉｄｔは、振動子ｉの信号成分について、損失に起因する複素振幅変化を表わす。

γ_ｉｎｊは、相互結合率を表わす。よって、（γ_ｉｎｊ／２）Ｊ_ｉｊＡ_ｊｄｔは、振動子ｉの信号成分について、振動子ｊから振動子ｉへの相互結合に起因する複素振幅変化を表わす。そして、（γ_ｉｎｊ／２）ΣＪ_ｉｊＡ_ｊｄｔは、振動子ｉの信号成分について、振動子ｉ以外のすべての振動子から振動子ｉへの相互結合に起因する複素振幅変化を表わす。

Ｄは、振動子の振幅の拡散係数を表わす。よって、√（Ｄ）ｄＷ_ｉは、振動子ｉの信号成分について、振幅拡散に起因する複素振幅変化を表わす。

複数の振動子のポテンシャル関数は、数１１により表わされる。ここで、Ｈ^〜（Ａ）は、ポテンシャル関数を表わし、太字のＡは、複数の振動子の信号成分について、複素振幅のベクトル表示を表わし、Ａ_ｉ ^＊は、Ａ_ｉの複素共役を表わす。

複数の振動子の時間発展を表わす複素ランジュバン方程式は、数１２のように書き換えられる。複数の振動子の信号成分について、複素振幅の定常状態の分布関数は、ポテンシャル関数を用いて、数１３により表わされ、より具体的には、数１４により表わされる。ここで、Ｐ_ｓｔ（Ａ）は、定常状態の分布関数を表わし、Ｃは、規格化因子を表わす。

ここで、複数の振動子の間の相互結合の強度が低いため、すべての振動子の強度が独立して数１５により表わされる定常状態の同じ強度に安定化する、と仮定する。なお、ｎ_ｓは、定常状態の同じ強度を表わす。すると、数１４の第１のｅｘｐの項のＤが小さいため、数１４の第１のｅｘｐの項が数１６により表わされるΠδ（｜Ａ_ｉ｜^２−ｎ_ｓ）に近似される。なお、δは、δ関数を表わし、Ｃ’は、新たな規格化因子を表わす。

ここで、振動子ｉの信号成分の複素振幅を、振幅√（ｎ_ｓ）及び位相θ_ｉを用いて、数１７のように表わす。すると、定常状態の分布関数は、数１８により表わされる。

さて、複数の振動子の系と複数のスピンの系を対応付けるために、数１９で表わされる疑似的な逆温度β＝１／（ｋ_ＢＴ）を定義する。ここで、ｋ_Ｂは、ボルツマン定数を表わし、Ｔは、複数のスピンの系の疑似的な温度を表わし、Ｄ_θ＝Ｄ／ｎ_ｓは、振動子の位相の拡散係数を表わす。

すると、定常状態の位相の分布関数は、数２０により表わされる。ここで、太字のθは、複数の振動子の信号成分について、位相のベクトル表示を表わし、Ｈ（θ）は、数２１により表わされる、ＸＹモデルのハミルトニアンを表わす。つまり、複数の振動子の系の定常状態の分布関数は、複数のスピンの系の疑似的な熱平衡状態のボルツマン分布に対応するのである。

（第３実施形態：計算方法）
次に、第３実施形態の結合振動子系の計算装置の構成を図１４に示す。結合振動子系の計算装置Ｃ５は、イテレーション部５１及びスピン系計算部５２から構成される。結合振動子系の計算装置Ｃ５に対して、結合振動子系の計算プログラムをインストールすることにより、イテレーション部５１及びスピン系計算部５２を実装することができる。

イテレーション部５１は、結合振動子系のダイナミクスを記述する微分方程式（数８）の各イテレーションにおいて、各振動子ｉ毎の微小時間後の振幅及び位相（数８のＡ_ｉの絶対値及び位相）を計算する。よって、イテレーション部５１は、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算することができる。

スピン系計算部５２は、イテレーション部５１により反復されたイテレーションの後に、微分方程式（数８）の定常状態における各振動子ｉ毎の位相（数８のＡ_ｉの位相）を計算することにより、結合振動子系の相互結合（数７のＪ_ｉｊ）に対応するＸＹモデル（数７）の平衡状態を計算する。よって、スピン系計算部５２は、ＸＹモデルの平衡状態を計算することにより、機械学習のボルツマンサンプリングを行なうことができる。

ここで、結合振動子系の全振動子数をＭとすると、数８の微分方程式のうち、右辺第１項の利得損失項及び右辺第３項の雑音項の計算量は、Ｍに比例するが、右辺第２項の相互結合項の計算量は、Ｍ^２に比例する。

しかし、数８の微分方程式のうち、右辺第２項の相互結合項の計算時間が、右辺第１項の利得損失項及び右辺第３項の雑音項の計算時間と比べて、同程度であるときには、図１４に示したように各イテレーションを行なえばよい。

つまり、ある回のイテレーションの後でのＡ_ｉを用いて計算した、右辺第１項の利得損失項、右辺第２項の相互結合項及び右辺第３項の雑音項を加算することにより、その回のイテレーションでの左辺第１項の微小時間変化を計算するのである。よって、右辺第２項の相互結合項を図３のような遅延（例えば、イテレーション数回分。）を伴わず計算しながら、右辺第１項の利得損失項及び右辺第３項の雑音項を高速に計算することができる。

ただし、数８の微分方程式のうち、右辺第２項の相互結合項の計算時間が、右辺第１項の利得損失項及び右辺第３項の雑音項の計算時間と比べて、十分に長いときには、図３と同じように各イテレーションを行なえばよい。

つまり、ある回のイテレーションの後でのＡ_ｉを用いて計算した右辺第１項の利得損失項及び右辺第３項の雑音項と、その回のイテレーションよりも数回（例えば、１回のみ。）遅れた回のイテレーションの後でのＡ_ｉを用いて計算した右辺第２項の相互結合項と、を加算することにより、その回のイテレーションでの左辺第１項の微小時間変化を計算するのである。よって、右辺第２項の相互結合項を余裕を持って計算しながら、右辺第１項の利得損失項及び右辺第３項の雑音項を高速に計算することができる。

結合振動子系の計算装置Ｃ５のうち、イテレーション部５１は、利得損失計算部５１１、相互結合計算部５１２、雑音計算部５１３及び振幅位相計算部５１４から構成される。

利得損失計算部５１１は、数８の微分方程式の各イテレーションにおいて、各振動子ｉ毎の利得及び損失に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、利得損失計算部５１１は、数８の微分方程式の各イテレーションにおいて、第ｎ回目のイテレーションの後でのＡ_ｉ ^（ｎ）を入力され、数８の右辺第１項を出力する。ここで、第１回目のイテレーションの前でのＡ_ｉ ^（０）は、０に初期設定される。

雑音計算部５１３は、数８の微分方程式の各イテレーションにおいて、各振動子ｉ毎の雑音に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、雑音計算部５１３は、数８の微分方程式の各イテレーションにおいて、第ｎ回目のイテレーションの後でのＡ_ｉ ^（ｎ）を入力され、数８の右辺第３項を出力する。

相互結合計算部５１２は、数８の微分方程式の各イテレーションにおいて、各振動子ｉの間の相互結合に起因する各振動子ｉ毎の振幅及び位相の微小時間変化を計算する。具体的には、相互結合計算部５１２は、数８の微分方程式の各イテレーションにおいて、第ｎ回目のイテレーションの後でのＡ_ｉ ^（ｎ）を入力され、数８の右辺第２項を出力する。ここで、相互結合計算部５１２での数８の右辺第２項の計算時間は、利得損失計算部５１１及び雑音計算部５１３での数８の右辺第１、３項の計算時間と比べて、同程度である。

振幅位相計算部５１４は、数８の微分方程式の各イテレーションにおいて、利得損失計算部５１１、雑音計算部５１３及び相互結合計算部５１２により各々計算された各振動子ｉ毎の振幅及び位相の微小時間変化に基づいて、各振動子ｉ毎の微小時間後の振幅及び位相を計算する。具体的には、振幅位相計算部５１４は、数８の微分方程式の各イテレーションにおいて、利得損失計算部５１１、雑音計算部５１３及び相互結合計算部５１２からの数８の右辺第１、３、２項を入力され、利得損失計算部５１１、雑音計算部５１３及び相互結合計算部５１２からの数８の右辺第１、３、２項を加算し、第（ｎ＋１）回目のイテレーションの後でのＡ_ｉ ^{（ｎ＋１）}を出力する。ここで、利得損失計算部５１１、雑音計算部５１３及び相互結合計算部５１２からの数８の右辺第１、３、２項は、第ｎ回目のイテレーションの後でのＡ_ｉ ^（ｎ）を用いて計算されたものである。

スピン系計算部５２は、イテレーション部５１により反復されたイテレーションの後に、数８の微分方程式の定常状態における各振動子ｉ毎の位相を計算することにより、結合振動子系の相互結合に対応するＸＹモデルの平衡状態を計算する。

このように、機械学習のボルツマンサンプリングを行なうために、機械学習を連続値スピンモデルへマッピングし、連続値スピンモデルの平衡状態を計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。そして、完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算するにあたり、電子回路で実装することにより、光学回路を安定化する必要をなくすことができる。

第３実施形態においても、第１実施形態の図６〜８と同様に、並列計算を行なうことができ、第１実施形態の図９〜１２と同様に、ビット量子化を行なうことができる。

本開示の結合振動子系の計算装置、プログラム及び方法は、次の目的に利用可能である：（１）組み合わせ最適化問題の最適解を求めるために、組み合わせ最適化問題をイジングモデルへマッピングし、イジングモデルの基底状態を計算するにあたり、光学回路を安定化する必要をなくす、（２）機械学習のボルツマンサンプリングを行なうために、機械学習を連続値スピンモデルへマッピングし、連続値スピンモデルの平衡状態を計算するにあたり、光学回路を安定化する必要をなくす、（３）完全結合系のニューロチップ等において、結合振動子系のダイナミクスを計算するにあたり、光学回路を安定化する必要をなくす。

Ｃ１、Ｃ２：イジングモデルの計算装置
Ｃ３、Ｃ４、Ｃ５：結合振動子系の計算装置
Ｐ１、Ｐ２、Ｐ３、Ｐ４：レーザパルス
１１、２１：レーザ発振部
１２、２２：リング共振器
１３、１５、１７、２４：相互結合実装部
１４、１６、１８：パルス遅延線
１９、２３：スピン測定部
３１、４１、５１：イテレーション部
３２、４２、５２：スピン系計算部
３３：ＦＩＦＯ
３１１、４１１、５１１：利得損失計算部
３１２、４１３、５１３：雑音計算部
３１３、４１２、５１２：相互結合計算部
３１４、４１４、５１４：振幅位相計算部

Claims (13)

1. 結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーション部と、
   前記イテレーション部により反復されたイテレーションの後に、前記微分方程式の定常状態における各振動子毎の位相を計算することにより、前記結合振動子系の相互結合に対応するスピン系の基底状態又は平衡状態を計算するスピン系計算部と、
   を備え、
   前記イテレーション部は、
   前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算部と、
   前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算部と、
   前記微分方程式の各イテレーションにおいて、前記利得損失計算部及び前記相互結合計算部により各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算部と、
   を備えることを特徴とする結合振動子系の計算装置。
2. 前記微分方程式の各イテレーションにおいて、各振動子毎の雑音に起因する各振動子毎の振幅及び位相の微小時間変化を計算する雑音計算部、
   をさらに備え、
   前記振幅位相計算部は、前記微分方程式の各イテレーションにおいて、前記利得損失計算部、前記相互結合計算部及び前記雑音計算部により各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する
   ことを特徴とする、請求項１に記載の結合振動子系の計算装置。
3. 結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーション部、
   を備え、
   前記イテレーション部は、
   前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算部と、
   前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算部と、
   前記微分方程式の各イテレーションにおいて、各振動子毎の雑音に起因する各振動子毎の振幅及び位相の微小時間変化を計算する雑音計算部と、
   前記微分方程式の各イテレーションにおいて、前記利得損失計算部、前記相互結合計算部及び前記雑音計算部により各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算部と、
   を備えることを特徴とする結合振動子系の計算装置。
4. 前記雑音計算部は、前記微分方程式の初期状態のみにおいて、各振動子毎の雑音に起因する各振動子毎の振幅及び位相の微小時間変化を計算する
   ことを特徴とする、請求項２又は３に記載の結合振動子系の計算装置。
5. 前記イテレーション部は、前回のイテレーションにおける各振動子毎の微小時間後の振幅及び位相が、前記微分方程式の定常状態に所定の精度で接近したときに、次回のイテレーションにおける各振動子毎の微小時間後の振幅及び位相の計算を停止する
   ことを特徴とする、請求項１から４のいずれかに記載の結合振動子系の計算装置。
6. 前記イテレーション部は、前記微分方程式の各イテレーションにおいて、各振動子毎の並列処理により、各振動子毎の微小時間後の振幅及び位相を計算する
   ことを特徴とする、請求項１から５のいずれかに記載の結合振動子系の計算装置。
7. 前記イテレーション部は、複数のイテレーション部を備え、
   前記複数のイテレーション部の各々は、前記微分方程式の各イテレーションにおいて、前記結合振動子系のうちの各々に特定のグループに属する振動子の微小時間後の振幅及び位相のデータを、前記複数のイテレーション部のその他との間で共有する
   ことを特徴とする、請求項６に記載の結合振動子系の計算装置。
8. 前回のイテレーションにおける前記データの共有処理が、次回のイテレーションにおける計算処理を停止／中断させないように、前記複数のイテレーション部の各々は、前記結合振動子系のうちの各々に特定のグループに属する振動子を割り当てられる
   ことを特徴とする、請求項７に記載の結合振動子系の計算装置。
9. 前記振幅位相計算部は、前記微分方程式の各イテレーションにおいて、所定の精度で量子化したうえで、各振動子毎の微小時間後の振幅及び位相を計算する
   ことを特徴とする、請求項１から８のいずれかに記載の結合振動子系の計算装置。
10. 前記振幅位相計算部により計算される各振動子毎の微小時間後の振幅についての量子化カットオフは、前記振幅位相計算部により計算される各振動子毎の微小時間後の振幅についての量子化の前記所定の精度を満足するように設定されており、かつ、前記微分方程式の定常状態における各振動子毎の振幅と比べて大きく設定されている
    ことを特徴とする、請求項９に記載の結合振動子系の計算装置。
11. 前記相互結合計算部により計算される各振動子毎の振幅の微小時間変化についての量子化ビット数は、前記振幅位相計算部により計算される各振動子毎の微小時間後の振幅についての量子化ビット数と比べて、少なくともｌｏｇ_２Ｍ（ただし、Ｍは、前記結合振動子系の全振動子数である。）だけ大きく設定されている
    ことを特徴とする、請求項９又は１０に記載の結合振動子系の計算装置。
12. 結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーションステップと、
    前記イテレーションステップにより反復されたイテレーションの後に、前記微分方程式の定常状態における各振動子毎の位相を計算することにより、前記結合振動子系の相互結合に対応するスピン系の基底状態又は平衡状態を計算するスピン系計算ステップと、
    をコンピュータに実行させ、
    前記イテレーションステップは、
    前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算ステップと、
    前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算ステップと、
    前記微分方程式の各イテレーションにおいて、前記利得損失計算ステップ及び前記相互結合計算ステップにより各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算ステップと、
    を備えることを特徴とする結合振動子系の計算プログラム。
13. 結合振動子系のダイナミクスを記述する微分方程式の各イテレーションにおいて、各振動子毎の微小時間後の振幅及び位相を計算するイテレーションステップと、
    前記イテレーションステップにより反復されたイテレーションの後に、前記微分方程式の定常状態における各振動子毎の位相を計算することにより、前記結合振動子系の相互結合に対応するスピン系の基底状態又は平衡状態を計算するスピン系計算ステップと、
    を備え、
    前記イテレーションステップは、
    前記微分方程式の各イテレーションにおいて、各振動子毎の利得及び損失に起因する各振動子毎の振幅及び位相の微小時間変化を計算する利得損失計算ステップと、
    前記微分方程式の各イテレーションにおいて、各振動子の間の相互結合に起因する各振動子毎の振幅及び位相の微小時間変化を計算する相互結合計算ステップと、
    前記微分方程式の各イテレーションにおいて、前記利得損失計算ステップ及び前記相互結合計算ステップにより各々計算された各振動子毎の振幅及び位相の微小時間変化に基づいて、各振動子毎の微小時間後の振幅及び位相を計算する振幅位相計算ステップと、
    を備えることを特徴とする結合振動子系の計算方法。

Images