特許 7582494 推定装置、推定方法、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-11-05

(45)【発行日】2024-11-13

(54)【発明の名称】推定装置、推定方法、及びプログラム

(51)【国際特許分類】

   G06N 99/00 20190101AFI20241106BHJP

【ＦＩ】

G06N99/00 180

【請求項の数】 7

(21)【出願番号】P 2023546600

(86)(22)【出願日】2021-09-07

(86)【国際出願番号】 JP2021032888

(87)【国際公開番号】W WO2023037417

(87)【国際公開日】2023-03-16

【審査請求日】2024-03-01

(73)【特許権者】

【識別番号】000004226

【氏名又は名称】日本電信電話株式会社

(74)【代理人】

【識別番号】110004381

【氏名又は名称】弁理士法人ＩＴＯＨ

(74)【代理人】

【識別番号】100107766

【弁理士】

【氏名又は名称】伊東 忠重

(74)【代理人】

【識別番号】100070150

【弁理士】

【氏名又は名称】伊東 忠彦

(74)【代理人】

【識別番号】100124844

【弁理士】

【氏名又は名称】石原 隆治

(72)【発明者】

【氏名】橋本 悠香

【審査官】山本 俊介

(56)【参考文献】

【文献】平松尚人 ほか，クープマンモード分解による部屋内温度勾配推定の実測データへの適用，第６２回システム制御情報学会 研究発表講演会 講演論文集 [CD-ROM]，一般社団法人 システム制御情報学会，2018年05月16日，pp.1-3（147-8）

【文献】石川 勲 ほか，ＲＫＨＳ上のＰｅｒｒｏｎ－Ｆｒｏｂｅｎｉｕｓ作用素を用いた力学系間の比較について，電子情報通信学会技術研究報告，日本，一般社団法人電子情報通信学会，2018年10月29日，第118巻，第284号，pp.175-182，ISSN：2432-6380

【文献】HASHIMOTO, Yuka ほか，Krylov Subspace Method for Nonlinear Dynamical Systems with Random Noise，Journal of Machine Learning Research [online]，vol. 21，2020年08月20日，pp.1-29（19-993），[検索日 2021.11.12], インターネット：<URL：https://www.jmlr.org/papers/volume21/19-993/19-993.pdf>

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０６Ｎ ３／００－９９／００

(57)【特許請求の範囲】

【請求項1】

複数の要素で構成される時系列データを入力として、前記時系列データからクープマン作用素を推定する作用素推定部と、
前記クープマン作用素を用いて、前記複数の要素の集団的な振動と前記要素間の相互作用とを表す位相モデルを推定する位相モデル推定部と、
を有する推定装置。

【請求項2】

前記位相モデル推定部は、
前記クープマン作用素を用いて、第１の最適化問題を勾配法により解き、
前記第１の最適化問題の解と前記クープマン作用素とを用いて、予め決められた所定の回数、第２の最適化問題を帰納的に解き、
前記第１の最適化問題の解と前記第２の最適化問題の解とを用いて、前記位相モデルを推定する、請求項１に記載の推定装置。

【請求項3】

前記時系列データをＸ（ｔ）＝［Ｘ_１（ｔ），・・・，Ｘ_Ｎ（ｔ）］、前記クープマン作用素をＫ、Ｂ_ｉ，ｋｕ＝ｕ_ｋｅ_ｉにより定義されるヒルベルト空間上の線形作用素をＢ_ｉ，ｋ、ｕ_ｋをベクトル値関数ｕの第ｋ成分、ｅ_ｉはｉ番目の要素のみ１でその他の要素は０であるＮ次元ベクトル、あるｔをｔ_０として、
前記第１の最適化問題は、

【数20】

と表される、請求項２に記載の推定装置。

【請求項4】

前記時系列データのデータ間隔をΔｔ、ある周波数をω∈［０，２π］、前記第１の最適化問題の解をλ＝ｅ＾（（√（－１））Δｔω）、ａ_ｉ，ｋ ^１、ｕ^１、前記クープマン作用素Ｋの固有値をλ_ｊ，ｉ＝ｅ＾（（√（－１））Δｔｊω）（ただし、ｊ＝２，・・・，Ｍ、ｉ＝１，・・・，Ｎ、Ｍは予め決められた２以上の整数）として、
前記第２の最適化問題は、

【数21】

と表され、
前記位相モデル推定部は、
ｊ＝２，・・・，Ｍに関して前記第２の最適化問題を帰納的に解く、請求項３に記載の推定装置。

【請求項5】

前記位相モデル推定部は、
前記第１の最適化問題の解λ、ａ_ｉ，ｋ ^１、ｕ^１と、前記第２の最適化問題の解ａ_ｉ，ｋ ^２，・・・，ａ_ｉ，ｋ ^Ｍとを用いて、位相結合関数を近似することで、前記周波数ωと前記位相結合関数とで構成される前記位相モデルを推定する、請求項４に記載の推定装置。

【請求項6】

複数の要素で構成される時系列データを入力として、前記時系列データからクープマン作用素を推定する作用素推定手順と、
前記クープマン作用素を用いて、前記複数の要素の集団的な振動と前記要素間の相互作用とを表す位相モデルを推定する位相モデル推定手順と、
をコンピュータが実行する推定方法。

【請求項7】

コンピュータを、請求項１乃至５の何れか一項に記載の推定装置として機能させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、推定装置、推定方法、及びプログラムに関する。

【背景技術】

【0002】

複数の要素で構成されるデータに関して、要素間の相互作用を調べることは、統計学、機械学習、物理学、分子動力学等といった様々な分野で共通に存在する課題である。物理学や分子動力学等においては、位相モデルと呼ばれるモデルにより、複数要素の集団的な振動とその間の相互作用の情報を記述する方法が提案されている（非特許文献１）。また、フーリエ級数やヒルベルト変換を用いて、与えられたデータから位相モデルを推定する方法が提案されている（非特許文献２）。位相モデルにおいては、位相関数と呼ばれる関数によりデータを位相の情報に変換した後、位相結合関数と呼ばれる関数により位相間の関係性を記述する。

【0003】

ここで、要素が１つの場合、位相モデルは１つの要素の振動のみを記述するが、この場合に関してはクープマン作用素を用いて位相モデルを推定する方法が提案されている（非特許文献３）。クープマン作用素は時系列データの時間発展を記述する線形作用素であり、ＲＫＨＳ（Reproducing Kernel Hilbert Space）やその一般化であるｖｖＲＫＨＳ（Vector-Valued RKHS、非特許文献４）を用いて、与えられた時系列データから、そのデータが従うクープマン作用素を推定する方法が提案されている（非特許文献５）。したがって、これらの方法を用いることで、要素数が１つの場合には、与えられたデータから位相モデルを推定することができる。

【先行技術文献】

【非特許文献】

【0004】

【文献】H. Nakao, S. Yasui, M. Ota, K. Arai, and Y. Kawamura, "Phase reduction and synchronization of a network of coupled dynamical elements exhibiting collective oscillations," Chaos 28, 045103 (2018).

【文献】B. Kralemann, et al., "In vivo cardiac phase response curve elucidates human respiratory heart rate variability," Nat. Commun. 4, 2418 (2013).

【文献】S. Shirasaka, W. Kurebayashi, H. Nakao, "Phase-amplitude reduction of transient dynamics far from attractors for limit-cycling systems," Chaos 27, 023119 (2017).

【文献】H. Q. Minh, L. Bazzani, and V. Murino, "A unifying framework in vector-valued reproducing kernel Hilbert spaces for manifold regularization and co-regularized multi-view learning," JMLR 17(25), 1-72 (2016).

【文献】Y. Hashimoto, I. Ishikawa, M. Ikeda, Y. Matsuo, and Y. Kawahara, "Krylov subspace method for nonlinear dynamical systems with random noise," JMLR 21(172), 1-29 (2020).

【発明の概要】

【発明が解決しようとする課題】

【0005】

フーリエ級数やヒルベルト変換を用いて位相モデルを推定する方法では位相関数全体を推定することは困難である一方で、クープマン作用素を用いて位相モデルを推定する方法では位相関数全体を推定することが可能である。しかしながら、クープマン作用素を用いて位相モデルを推定する方法は、要素が１つのみであることが必要である。

【0006】

本発明の一実施形態は、上記の点に鑑みてなされたもので、クープマン作用素を用いて複数の要素に対する位相モデルを推定することを目的とする。

【課題を解決するための手段】

【0007】

上記目的を達成するため、一実施形態に係る推定装置は、複数の要素で構成される時系列データを入力として、前記時系列データからクープマン作用素を推定する作用素推定部と、前記クープマン作用素を用いて、前記複数の要素の集団的な振動と前記要素間の相互作用とを表す位相モデルを推定する位相モデル推定部と、を有する。

【発明の効果】

【0008】

クープマン作用素を用いて複数の要素に対する位相モデルを推定することができる。

【図面の簡単な説明】

【0009】

【図1】本実施形態に係る推定装置のハードウェア構成の一例を示す図である。

【図2】本実施形態に係る推定装置の機能構成の一例を示す図である。

【図3】本実施形態に係る位相モデル推定処理の一例を示すフローチャートである。

【図4】クープマン作用素の固有値を複素平面上にプロットした散布図の一例を示す図（その１）である。

【図5】観測データから推定した位相関数の一例を示す図である。

【図6】ＦＨＮモデルから計算した変換の一例を示す図である。

【図7】相互作用の反対称部の値の一例を示す図である。

【図8】クープマン作用素の固有値を複素平面上にプロットした散布図の一例を示す図（その２）である。

【図9】相互作用の強さを表すヒートマップの一例を示す図である。

【発明を実施するための形態】

【0010】

以下、本発明の一実施形態について説明する。本実施形態では、クープマン作用素を用いて複数の要素に対する位相モデルを推定することができる推定装置１０について説明する。

【0011】

＜理論的構成＞
まず、本実施形態に係る推定装置１０がクープマン作用素を用いて複数の要素に対する位相モデルを推定する際の理論的構成について説明する。

【0012】

≪１．設定≫
以下のように生成された時系列データを解析することを考える。

【0013】

【数1】

とする。なお、ｄは予め決められた自然数である。以下、明細書のテキスト中では、上記の数１に示すＸ（カリグラフィー文字のＸ）を「Ｘｓｅｔ」と表記する。また、他の記号と混同が生じない場合には、明細書のテキスト中では黒板太字（中抜き文字）を通常の文字で表記する（例えば、上記の数１のｄ次元実数空間はＲ^ｄと表記する。）。

【0014】

以下の式（１）に示すネットワーク（結合力学系のネットワークモデル）で記述されるＮ個の要素をＸ_１，・・・，Ｘ_Ｎ∈Ｘｓｅｔとする。

【0015】

【数2】

ただし、初期値はＸ_ｉ（０）＝ｘ_ｉ，０（ｉ＝１，・・・，Ｎ）である。Ｆ_ｉは要素Ｘ_ｉの個別のダイナミクスを表し、Ｇ_ｉ，ｋは要素Ｘ_ｉに対する要素Ｘ_ｋの影響を表す。

【0016】

このとき、Ｘ＝［Ｘ_１，・・・，Ｘ_Ｎ］とおく。上記の式（１）に示すモデルによって、一定の時間間隔Δｔでデータが生成される、つまり、ｘ_ｉ，ｌ＝Ｘ_ｉ（Δｔ・ｌ）という時系列データが生成されるとする。これらの時系列データが観測可能なデータであり、観測データとも呼ぶ。なお、ｌは小文字のＬである。

【0017】

Ｎ個の要素を支配している共通の周波数ωが存在し、更に、各要素は他の要素と弱く相互作用しているものと仮定する。つまり、上記の式（１）に示すモデルは、以下の式（２）に示す位相モデルに縮約されると仮定する。

【0018】

【数3】

ここで、θ_ｉ∈［０，２π）はＸ_ｉに対する位相変数、Γ_ｉは要素Ｘ_ｉに影響を与える要素間の相互作用を表す。Ｘ_ｉをθ_ｉ∈［０，２π）に変換する関数を位相関数と呼び、Γ_ｉを位相結合関数とも呼ぶ。本実施形態では、観測データのみから（つまり、Ｆ_ｉとＧ_ｉ，ｋを知らない状態で）、位相関数、周波数ω、相互作用Γ_ｉを推定することを考える。

【0019】

≪２．クープマン作用素を用いた位相モデルの推定≫
Ｘｓｅｔ^ＮからＣ^Ｎ（Ｎ次元複素数空間）への関数ｖであって、ｖ（ｘ）∈Ｃ^Ｎの第ｉ成分がｘ∈Ｘｓｅｔ^Ｎの第ｉ成分のみに依存する関数ｖで構成されるヒルベルト空間をＨとする。例えば、以下の行列値カーネル関数Φ：Ｘｓｅｔ^Ｎ×Ｘｓｅｔ^Ｎ→Ｃ^Ｎ×Ｎから生成されるｖｖＲＫＨＳをＨとして設定することが可能である。

【0020】

［Φ（ｘ_１，ｘ_２）］_ｉ，ｊ＝ｋ（（ｘ_１，ｉ，ｉ），（ｘ_２，ｊ，ｊ））
ただし、ｋは

【0021】

【数4】

上の複素数値正定値カーネル、ｘ_１，ｉはｘ_１∈Ｘｓｅｔ^Ｎの第ｉ成分、ｘ_２，ｊはｘ_２∈Ｘｓｅｔ^Ｎの第ｊ成分を表す。なお、ｖｖＲＫＨＳの構成の詳細については、例えば、非特許文献４を参照されたい。

【0022】

ｖ∈Ｈに対して、Ｈ上のクープマン作用素ＫをＫｖ（Ｘ（ｔ））＝ｖ（Ｘ（ｔ＋Δｔ））を満たすような線形作用素として定義する。ある自然数ｌ（ｌは小文字のＬ）に対して、ｔ_０＝Δｔ・ｌとする。このとき、以下の式（３）に示す最適化問題を考える。

【0023】

【数5】

ただし、Ｂ_ｉ，ｋはＢ_ｉ，ｋｕ＝ｕ_ｋｅ_ｉにより定義されるＨ上の線形作用素で、ｕ_ｋはベクトル値関数ｕの第ｋ成分、ｅ_ｉはｉ番目の要素のみ１でその他の要素は０であるＮ次元ベクトルである。なお、この式（３）に示す最適化問題は、例えば、勾配法等により解くことが可能である。

【0024】

あるω∈［０，２π］により、

【0025】

【数6】

と表せるものとし、更に、

【0026】

【数7】

に値が近いクープマン作用素Ｋの固有値が各ｊに対してＮ個存在するものとする。以下、これらの固有値をλ_ｊ，ｉ（ｉ＝１，・・・，Ｎ）と表す。なお、Ｍは予め決められた２以上の自然数である。

【0027】

固有値λ_ｊ，ｉに対する固有ベクトルが存在するものとし、これをｖ_ｊ，ｉと表す。このとき、各ｊ＝２，・・・，Ｍに対してｖ_ｊ，１（Ｘ（ｔ_０）），・・・，ｖ_ｊ，Ｎ（Ｘ（ｔ_０））が一次独立であるものとすると、

【0028】

【数8】

を満たすｃ_ｊ，ｉ∈Ｃが存在する。ただし、上記の数８の右辺の１はすべての成分が１のＮ次元ベクトルである。

【0029】

このとき、ｊ＝２，・・・，Ｍに対して、以下の式（４）に示す最適化問題Ｐ_ｊを帰納的に考える。

【0030】

【数9】

ただし、λ、ａ_ｉ，ｋ ^１、ｕ^１は上記の式（３）に示す最適化問題の解、ｌ＜ｊに対してａ_ｉ，ｋ ^ｌは最小化問題Ｐ_ｌの解、ｌ＜ｊに対して

【0031】

【数10】

である。なお、ｌは小文字のＬである。また、

【0032】

【数11】

は成分ごとの積を表す。

【0033】

上記の式（４）に示す最適化問題Ｐ_ｊを帰納的にｊ＝Ｍまで解くと、以下が得られる。なお、上記の式（４）に示す最適化問題は線形な問題であるため、解析的に解くことが可能である。

【0034】

【数12】

ただし、ｕ_ｉ，ｊはｕ_ｊのｉ番目の要素である。

【0035】

θ_ｉ，ｊ（ｔ）＝ａｒｇ（ｕ_ｉ，ｊ（Ｘ_ｉ（ｔ）））、ｒ_ｉ，ｊ（ｔ）＝｜ｕ_ｉ，ｊ（Ｘ_ｉ（ｔ））｜とおく。ｕ_ｊの定義と式（４）より、

【0036】

【数13】

である。よって、以下が成立する。

【0037】

【数14】

また、ｖ_ｉ，ｊはクープマン作用素Ｋの固有値

【0038】

【数15】

に対する固有ベクトルであるため、ｕ_ｊの定義より、

【0039】

【数16】

となる。よって、

【0040】

【数17】

となるため、θ_ｉ，ｊ（ｔ）をｊθ_ｉ，１（ｔ）で近似すると、以下の式（５）が成立する。

【0041】

【数18】

したがって、上記の式（５）において、

【0042】

【数19】

とおけば上記の式（２）が得られる。なお、θ_ｉ，１（ｔ）＝ａｒｇ（ｕ_ｉ，１（Ｘ_ｉ（ｔ）））が位相関数である。

【0043】

≪３．クープマン作用素の推定≫
もし観測データからクープマン作用素Ｋを推定することができれば、「２．クープマン作用素を用いた位相モデルの推定」で説明した方法により、当該観測データが従う位相モデルを推定することが可能である。観測データからクープマン作用素Ｋを推定する方法としては、例えば、非特許文献５に記載されている方法を用いればよい。

【0044】

＜推定装置１０のハードウェア構成＞
次に、本実施形態に係る推定装置１０のハードウェア構成について、図１を参照しながら説明する。図１に示すように、本実施形態に係る推定装置１０は一般的なコンピュータ又はコンピュータシステムのハードウェア構成で実現され、入力装置１０１と、表示装置１０２と、外部Ｉ／Ｆ１０３と、通信Ｉ／Ｆ１０４と、プロセッサ１０５と、メモリ装置１０６とを有する。これらの各ハードウェアは、それぞれがバス１０７により通信可能に接続される。

【0045】

入力装置１０１は、例えば、キーボード、マウス、タッチパネル、各種物理ボタン等である。表示装置１０２は、例えば、ディスプレイ、表示パネル等である。なお、推定装置１０は、例えば、入力装置１０１及び表示装置１０２のうちの少なくとも一方を有していなくてもよい。

【0046】

外部Ｉ／Ｆ１０３は、記録媒体１０３ａ等の外部装置とのインタフェースである。推定装置１０は、外部Ｉ／Ｆ１０３を介して、記録媒体１０３ａの読み取りや書き込み等を行うことができる。なお、記録媒体１０３ａとしては、例えば、ＣＤ（Compact Disc）、ＤＶＤ（Digital Versatile Disk）、ＳＤメモリカード（Secure Digital memory card）、ＵＳＢ（Universal Serial Bus）メモリカード等が挙げられる。

【0047】

通信Ｉ／Ｆ１０４は、推定装置１０を通信ネットワークに接続するためのインタフェースである。プロセッサ１０５は、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）等の各種演算装置である。メモリ装置１０６は、例えば、ＨＤＤ（Hard Disk Drive）、ＳＳＤ（Solid State Drive）、フラッシュメモリ、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）等の各種記憶装置である。

【0048】

本実施形態に係る推定装置１０は、図１に示すハードウェア構成を有することにより、後述する各種処理を実現することができる。なお、図１に示すハードウェア構成は一例であって、推定装置１０は、例えば、複数のプロセッサ１０５を有していてもよいし、複数のメモリ装置１０６を有していてもよいし、図示したハードウェア以外の様々なハードウェアを有していてもよい。

【0049】

＜推定装置１０の機能構成＞
次に、本実施形態に係る推定装置１０の機能構成について、図２を参照しながら説明する。図２に示すように、本実施形態に係る推定装置１０は、クープマン作用素推定部２０１と、位相モデル推定部２０２と、記憶部２０３とを有する。クープマン作用素推定部２０１及び位相モデル推定部２０２は、例えば、推定装置１０にインストールされた１以上のプログラムが、プロセッサ１０５に実行させる処理により実現される。また、記憶部２０３は、例えば、メモリ装置１０６により実現される。なお、記憶部２０３は、推定装置１０と通信ネットワークを介して接続される記憶装置（例えば、データベースサーバ等）により実現されてもよい。

【0050】

クープマン作用素推定部２０１は、観測データからクープマン作用素Ｋを推定する。ここで、クープマン作用素推定部２０１には、データ取得部２１１と、作用素推定部２１２とが含まれている。データ取得部２１１は、記憶部２０３から観測データを取得する。作用素推定部２１２は、例えば、非特許文献５に記載されている方法を用いて、データ取得部２１１によって取得された観測データからクープマン作用素Ｋを推定する。

【0051】

位相モデル推定部２０２は、クープマン作用素推定部２０１によって推定されたクープマン作用素Ｋを用いて、「２．クープマン作用素を用いた位相モデルの推定」で説明した方法により、当該観測データが従う位相モデルを推定する。

【0052】

記憶部２０３は、観測データが記憶されている。なお、記憶部２０３には、位相モデル推定部２０２によって推定された位相モデル（位相関数、周波数ω、相互作用Γ_ｉ）が記憶されてもよい。

【0053】

＜位相モデル推定処理＞
次に、本実施形態に係る位相モデル推定処理について、図３を参照しながら説明する。

【0054】

クープマン作用素推定部２０１のデータ取得部２１１は、記憶部２０３から観測データを取得する（ステップＳ１０１）。

【0055】

クープマン作用素推定部２０１の作用素推定部２１２は、上記のステップＳ１０１で取得された観測データからクープマン作用素Ｋを推定する（ステップＳ１０２）。なお、上述したように、作用素推定部２１２は、例えば、非特許文献５に記載されている方法を用いて、観測データからクープマン作用素Ｋを推定すればよい。

【0056】

位相モデル推定部２０２は、上記のステップＳ１０２で推定されたクープマン作用素Ｋを用いて、上記の式（３）に示す最適化問題を解く（ステップＳ１０３）。なお、位相モデル推定部２０２は、例えば、勾配法等により上記の式（３）に示す最適化問題を解くことができる。

【0057】

位相モデル推定部２０２は、ｊ←２とする（ステップＳ１０４）。

【0058】

位相モデル推定部２０２は、ｊ＜Ｍ＋１であるか否かを判定する（ステップＳ１０５）。なお、Ｍは予め決められた２以上の自然数である。

【0059】

上記のステップＳ１０５でｊ＜Ｍ＋１であると判定された場合、位相モデル推定部２０２は、上記のステップＳ１０２で推定されたクープマン作用素Ｋを用いて、上記の式（４）に示す最適化問題を解く（ステップＳ１０６）。なお、上記の式（４）に示す最適化問題は線形な問題であるため、位相モデル推定部２０２は、この式（４）に示す最適化問題を解析的に解くことができる。

【0060】

位相モデル推定部２０２は、ｊ←ｊ＋１として（ステップＳ１０７）、上記のステップＳ１０５に戻る。これにより、ｊ＝２，・・・，Ｍに対して、上記の式（４）に示す最適化問題が再帰的に解かれることになる。

【0061】

上記のステップＳ１０５でｊ＜Ｍ＋１であると判定されなかった場合、位相モデル推定部２０２は、上記の式（３）に示す最適化問題の解と上記の式（４）に示す最適化問題の解とを用いて、上記の式（５）により位相モデルを推定する（ステップＳ１０８）。これにより、上記のステップＳ１０１で取得された観測データが従う位相モデルが推定されたことになる。

【0062】

＜評価＞
以下、本実施形態に係る推定装置１０によって推定された位相モデルの評価について説明する。

【0063】

≪位相関数、位相結合関数の推定≫
Ｘｓｅｔ＝Ｒ^２上のＦＨＮ（FitzHugh-Nagumo）モデルから初期値を変化させることにより１０個の時系列データ（観測データ）ｘ_０ ^ｉ，・・・，ｘ_２０００ ^ｉ（ｉ＝１，・・・，１０）を生成し、クープマン作用素Ｋを推定した上で、「２．クープマン作用素を用いた位相モデルの推定」で説明した方法により、Ｘ_ｉ（ｔ）からθ_ｉ（ｔ）への変換を表す位相関数η_ｉ（η_ｉ（ｘ）＝ａｒｇ（ｕ_ｉ，１（ｘ）），ｘ∈Ｒ^２）と、上記の式（５）における相互作用Γ_ｉ（ψ_１，・・・，ψ_Ｎ）とを計算した。

【0064】

ＦＨＮモデルは、上記の式（１）に示すネットワークにおいて、Ｘ_ｉ（ｔ）＝［ｙ_ｉ（ｔ），ｚ_ｉ（ｔ）］∈Ｒ^２に対して、Ｆ_１（Ｘ_ｉ）＝Ｆ（Ｘ_ｉ）＝［ｙ_ｉ（ｙ_ｉ－ｃ）（１－ｙ_ｉ）－ｚ_ｉ，μ^－１（ｙ_ｉ－ｄｚ_ｉ）］、Ｇ_ｉ，ｊ（Ｘ_ｉ，Ｘ_ｊ）＝Ｇ（Ｘ_ｉ，Ｘ_ｊ）＝［０．０１（ｚ_ｉ－ｚ_ｊ），０］、ｃ＝－０．１、ｄ＝０．５、μ＝１００としたものである。

【0065】

このとき、クープマン作用素推定部２０１によって推定されたクープマン作用素Ｋの固有値を複素平面上にプロットした散布図を図４に示す。なお、横軸が実部、縦軸が虚部である。

【0066】

図４に示すように、値のほぼ等しい２つの固有値が重なっている箇所が３箇所（Ｑ_１１、Ｑ_１２、Ｑ_１３）ある。これら３箇所のうち１に最も近い箇所（Ｑ_１３）に存在する２つの固有値をλ_１，１、λ_１，２とする。また、固有値λ_１，１に対する固有ベクトルをｖ_１，１、固有値λ_１，２に対する固有ベクトルをｖ_１，２とする。このとき、ｓ＝１８９９として、ｃ_１，１Σ_ｉ＝１ ^１０ｖ_１，１（ｘ_ｓ ^ｉ）／１０＋ｃ_１，２Σ_ｉ＝１ ^１０ｖ_１，２（ｘ_ｓ ^ｉ）／１０＝１を満たすｃ_１，１，ｃ_１，２∈Ｃに対して、ｕ_１ ^０＝ｃ_１，１ｖ_１，１＋ｃ_１，２ｖ_１，２とおき、上記の式（３）に示す最適化問題をλ＝λ_１，１、ｕ＝ｕ_１ ^０、ａ_１，２＝ａ_２，１＝０．００２、ａ_１，１＝ａ_２，２＝０を初期値として解いた。更に、Ｍ＝３とし、値のほぼ等しい２つの固有値が重なっている箇所のうち、λ_１，１及びλ_１，２が存在する箇所以外の箇所に存在する固有値をそれぞれλ_２，ｉ（ｉ＝１，２）、λ_３，ｉ（ｉ＝１，２）と設定し、上記の式（４）に示す最適化問題を解いた。

【0067】

この例ではＮ＝２であり、上記の式（１）においてＦ_１＝Ｆ_２、Ｇ_１，２＝Ｇ_２，１であるため、η_１＝η_２となるはずである。そこで、位相モデル推定部２０２によってη_１のみを計算した。その結果を図５に示す。一方で、Ｘ_ｉ（ｔ）からθ_ｉ（ｔ）への変換をＦＨＮモデルから直接計算した結果を図６に示す。図５と図６を比較すると、観測データから推定した結果（図５）は、ＦＨＮモデルから直接計算した結果（図６）と近いことがわかる。

【0068】

また、相互作用Γ_ｉの反対称部Γ_α（ψ）＝Γ_１（ψ）－Γ_２（－ψ）を計算した結果を図７に示す。これに関しても、図７に示すように、観測データから推定した値と、ＦＨＮモデルから直接計算した値とが近いことがわかる。

【0069】

≪相互作用の強さの推定≫
Ｘｓｅｔ＝Ｒ^２上のＳＬモデルから初期値を変化させることにより１０個の時系列データ（観測データ）ｘ_０ ^ｉ，・・・，ｘ_２０００ ^ｉ（ｉ＝１，・・・，１０）を生成し、クープマン作用素Ｋを推定した上で、上記の式（３）に示す最適化問題を解き、相互作用の強さを表すａ_ｉ，ｋを推定した。

【0070】

ＳＬモデルは、上記の式（１）に示すネットワークにおいて、Ｘ_ｉ（ｔ）＝［ｙ_ｉ（ｔ），ｚ_ｉ（ｔ）］∈Ｒ^２に対して、Ｆ_１（Ｘ_ｉ）＝Ｆ（Ｘ_ｉ）＝［ｙ_ｉ－ａｚ_ｉ－（ｙ_ｉ ^２＋ｚ_ｉ ^２）（ｂｙ_ｉ＋ｚ_ｉ），ｚ_ｉ－（ｙ_ｉ ^２＋ｚ_ｉ ^２）（ｂｙ_ｉ＋ｚ_ｉ）］、Ｇ（Ｘ_ｉ，Ｘ_ｊ）＝［０．０１（ｚ_ｊ－ｚ_ｉ），０］、Ｇ_１，２（Ｘ_ｉ，Ｘ_ｊ）＝Ｇ_２，１（Ｘ_ｉ，Ｘ_ｊ）＝Ｇ（Ｘ_ｉ，Ｘ_ｊ）、Ｇ_１，３（Ｘ_ｉ，Ｘ_ｊ）＝Ｇ_２，３（Ｘ_ｉ，Ｘ_ｊ）＝Ｇ_３，１（Ｘ_ｉ，Ｘ_ｊ）＝Ｇ_３，２（Ｘ_ｉ，Ｘ_ｊ）＝０．１Ｇ（Ｘ_ｉ，Ｘ_ｊ）、ａ＝２、ｂ＝１としたものである。

【0071】

このとき、クープマン作用素推定部２０１によって推定されたクープマン作用素Ｋの固有値を複素平面上にプロットした散布図を図８に示す。なお、横軸が実部、縦軸が虚部である。

【0072】

図８に示すように、値のほぼ等しい３つの固有値が重なっている箇所が３箇所（Ｑ_２１、Ｑ_２２、Ｑ_２３）ある。これら３箇所Ｑ_２１、Ｑ_２２、Ｑ_２３に存在する固有値をそれぞれλ_１，λ_２，λ_３とする。また、固有値λ_１に対する固有ベクトルをｖ_１、固有値λ_２に対する固有ベクトルをｖ_２、固有値λ_３に対する固有ベクトルをｖ_３とする。このとき、ｓ＝１８９９として、ｃ_１Σ_ｉ＝１ ^１０ｖ_１（ｘ_ｓ ^ｉ）／１０＋ｃ_２Σ_ｉ＝１ ^１０ｖ_２（ｘ_ｓ ^ｉ）／１０＋ｃ_３Σ_ｉ＝１ ^１０ｖ_３（ｘ_ｓ ^ｉ）／１０＝１を満たすｃ_１，ｃ_２，ｃ_３∈Ｃに対して、ｕ_１ ^０＝ｃ_１ｖ_１＋ｃ_２ｖ_２＋ｃ_３ｖ_３とした。

【0073】

この例ではａ_ｉ，ｋの大きさを求めることのみが目的であるため、上記の式（３）において、λ＝λ_１、ｕ＝ｕ_１ ^０と固定し、線形な問題として式（３）に示す最適化問題を解いた。これにより求められたａ_ｉ，ｋの大きさを表すヒートマップを図９に示す。図９に示すように、要素Ｘ_１と要素Ｘ_２の間の相互作用が大きいと推定されており、これは、元のモデル（ＳＬモデル）で要素Ｘ_１と要素Ｘ_２の間の相互作用が大きいことを反映した結果となっている。

【0074】

＜まとめ＞
以上のように、本実施形態に係る推定装置１０は、クープマン作用素を用いて複数の要素に対する位相モデルを精度良く推定することができる。このため、本実施形態に係る推定装置１０によって推定された位相モデルを用いて、複数の要素で構成されるデータに関してそれらの要素間の相互作用を抽出する等といった解析を行うことが可能となる。なお、本実施形態に係る推定装置１０は、例えば、推定した位相モデルを用いて上記の解析を行ってもよいし、この解析結果を用いて結合力学系のネットワークモデルで表されるシステム（又は、そのシステムを構成する機器等）を制御してもよい。

【0075】

本発明は、具体的に開示された上記の実施形態に限定されるものではなく、請求の範囲の記載から逸脱することなく、種々の変形や変更、既知の技術との組み合わせ等が可能である。

【符号の説明】

【0076】

１０ 推定装置
１０１ 入力装置
１０２ 表示装置
１０３ 外部Ｉ／Ｆ
１０３ａ 記録媒体
１０４ 通信Ｉ／Ｆ
１０５ プロセッサ
１０６ メモリ装置
１０７ バス
２０１ クープマン作用素推定部
２０２ 位相モデル推定部
２０３ 記憶部
２１１ データ取得部
２１２ 作用素推定部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】