はじめに
Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Controlの第7.4節、クープマン作用素論について紹介する。普通に理解するのにてこずった。
概要
クープマン作用素論は、非線形ダイナミクスを線形作用素として解析する数学的フレームワークであり、有限次元の非線形システムを無限次元の線形系として扱うことを可能にする。この理論は、流体力学、データ駆動モデリング、機械学習など、多くの応用分野で活用されている。
一方、以前投稿した動的モード分解は、クープマン作用素を近似的に数値計算するデータドリブンな手法であり、観測データから固有値・固有関数・モードを抽出する。これは、複雑なシステムの支配的なダイナミクスを抽出する上で強力なツールとなる。
本稿では、クープマン作用素論を包括的に理解し、理論的基盤と数値解析の双方の観点からその役割を明らかにすることを目指す。
クープマン作用素論への入門
動的システム（Dynamical System）の導入
こちらの記事で用いた定義と同じものを用いる。先ずは、非線形常微分方程式で表現される連続力学系を考える。
時間 $t\in\mathbb{R}$ 、 $\boldsymbol{x}\in M$ は状態空間、 $\boldsymbol{f}$ は非線形のベクトル場（次の $\boldsymbol{F}$ の流れ場）として：
\cfrac{d\boldsymbol{x}}{dt}
=\boldsymbol{f}(\boldsymbol{x},t,\beta)
\tag{1}
これを元に離散力学系の基本定義は以下のように導かれる:
\begin{align}
\boldsymbol{x}(t)
&=\boldsymbol{F}\left(\boldsymbol{x}(t-1)\right)
\tag{2} \\
&=\boldsymbol{F}^t\left(\boldsymbol{x}(0)\right)
\tag{3}
\end{align}
参考書内では、 $\boldsymbol{x}_{k+1}=\boldsymbol{F}(\boldsymbol{x_k})$ のように記述されているが、Steven L. Brunton et al, 2022に倣って本稿では上のような記述を採用する。
クープマン理論の数学的記述
クープマン作用素は、非線形系の状態空間で定義される観測関数をダイナミクスと共に時間発展させる線形作用素であり、次式のように与えられる。
\begin{align}
&&\mathcal{K}^tg(\boldsymbol{x})
&=g\left(\boldsymbol{F}^t(\boldsymbol{x})\right)
\tag{4} \\
&\Rightarrow&
\mathcal{K}^tg
&=g\circ\boldsymbol{F}^t
\tag{5}
\end{align}
$g:M\rightarrow\mathbb{R}$ は「観測関数」と呼ばれる実測値観測関数であり、その関数空間を $\mathcal{G}(M)$ と記述する。
この $(5)$ 式は、クープマン作用素が時間 $t$ による観測関数の進化 $t\rightarrow g,\ g\in\mathcal{G}(M)$ と解釈する事ができる：
g_t:=\mathcal{K}^tg,\quad g_0:=g
\tag{6}
クープマン作用素 $\mathcal{K}$ は、状態空間 $M$ の各点 $\boldsymbol{x}$ に作用するのではなく、状態空上の関数 $g$ に作用する。 また、この $g$ の関数空間 $\mathcal{G}$ は次のような数学的性質を持つことが求められる。
-
ベクトル空間
- 関数$g_1,\ g_2\in\mathcal{G}*$に対し、加法性とスカラー倍が成立する。
- すなわち、$\alpha g_1+\beta g_2\in\mathcal{G}$が成り立つ。
-
内積が定義できる
-
完備性を持つ
-
$M$上の連続関数のような特定の関数を含む
この関数空間 $\mathcal{G}$ が線形（ベクトル）空間である場合、クープマン作用素も線形となる。
\begin{align}
\mathcal{K}^t
\left(
\alpha_1g_1(\boldsymbol{x})+
\alpha_2g_2(\boldsymbol{x})
\right)
&=
\alpha_1g_1
\left(
\boldsymbol{F}^t(\boldsymbol{x})
\right)
+
\alpha_2g_2
\left(
\boldsymbol{F}^t(\boldsymbol{x})
\right)
\tag{7} \\
&=
\alpha_1\mathcal{K}^tg_1(\boldsymbol{x})
+
\alpha_2\mathcal{K}^tg_2(\boldsymbol{x})
\tag{8}
\end{align}
このクープマン作用素の線形性は元のダイナミクス（流れ場） $\boldsymbol{F}^t:M\rightarrow M$ が非線形であっても成立する。
即ち、クープマンフレームワークは、有限次元の状態空間 $M$ を無限次元の関数空間 $\mathcal{G}(M)$ に拡張することで、非線形なシステムに対しても線形な表現を得られる、というものである。
離散時間におけるクープマン力学系
離散時間 $t\in\mathbb{N}$ の場合、時間発展 $\boldsymbol{F}^t$ は、写像 $\boldsymbol{F}\equiv\boldsymbol{F}^1$ を $t$ 回繰り返したものになる：
\boldsymbol{F}^t(\boldsymbol{x})
=\boldsymbol{F}(\boldsymbol{F}(\cdots(\boldsymbol{F}(\boldsymbol{x}))))
\tag{9}
クープマン作用素もこれと同様に、作用素 $\mathcal{K}\equiv\mathcal{K}^1$ を繰り返し適用することで得られる。このようにして得られる $\mathcal{K}$ は、(可算個の) 繰り返し合成によって生成される半群の生成作用素と解釈される。
結果として、力学系の更新：
\boldsymbol{x}_{k+1}=\boldsymbol{F}(\boldsymbol{x}_k)\tag{10}
に対応する、観測関数の更新を得る。：
g_{k+1}=\mathcal{K}g_k
\tag{11}
連続時間におけるクープマン力学系
連続時間系の場合、ダイナミクス（流れ場）は半群の性質を持つ：
\boldsymbol{F}^{t+s}
=\boldsymbol{F}^t(\boldsymbol{F}^s(\boldsymbol{x}))
\quad
\forall\boldsymbol{x},\ \lbrace t,s\rbrace\in\mathbb{R_{\geq0}}
\tag{12}
つまり、時間 $s$ 進めた後にさらに $t$ 進めるのは、直接 $t+s$ 進めるのと同じになる。もしここで、写像 $\boldsymbol{F}^t$ が逆射を持てば、これは群となる。
そして、クープマン作用素 $\mathcal{K}^t$ もこれと同じ性質を持つ：
\mathcal{K}^{t+s}=\mathcal{K}^t\mathcal{K}^s
\tag{13}
さらに、（ややこしいのだけれども）連続時間系ではクープマン作用素の時間微分（無限小生成作用素）を考えることができ、これをリー作用素 $\mathcal{L}$ として定義する：
\mathcal{L}g
:=\lim_{t\to0}\frac{\mathcal{K}^tg-g}{t}
=\lim_{t\to0}\frac{g\circ\boldsymbol{F}^t-g}{t}
\tag{14}
このリー作用素は連続力学系 $d\boldsymbol{x}/dt
=\boldsymbol{f}(\boldsymbol{x},t,\beta)$ の観測関数 $g(\boldsymbol{x})$ を時間微分することでも得られる：
\begin{align}
\cfrac{d}{dt}\,g(\boldsymbol{x})
&=\nabla\boldsymbol{g}\cdot\cfrac{d\boldsymbol{x}}{dt}
\tag{15} \\
&=\nabla\boldsymbol{g}\cdot\boldsymbol{f}(\boldsymbol{x},t,\beta)
\tag{16}
\end{align}
また、この左辺については以下のように変形できる：
\cfrac{d}{dt}\,g(\boldsymbol{x})
=\lim_{\tau\to0}\cfrac{g(\boldsymbol{x}(t+\tau))-g(\boldsymbol{x}(t))}{\tau}
=\mathcal{L}(g(\boldsymbol{x}(t)))
\tag{17}
故に
\mathcal{L}g=\nabla g\cdot\boldsymbol{f}(\boldsymbol{x},t,\beta)
\tag{18}
（余談）共役作用素について
ヒルベルト空間上の線形作用素 $A$ に対して、ある作用素 $A^$ が存在し、以下を満たす場合に、$A^$ を $A$ の共役作用素と呼ぶ：
\langle Ax, y\rangle=\langle x, A^*y\rangle
\tag{19}
良く知られた例で言えば、行列の場合の共役作用素はエルミート共役であり、以下の二つの行列は定義 $(19)$ を満足し、共役関係にある。
A=
\begin{bmatrix}
1 & i \\ 1 & 2
\end{bmatrix}
,\quad
A^*=
\begin{bmatrix}
1 & 1 \\ -i & 2
\end{bmatrix}
\tag{20}
リー作用素の共役作用素はリウヴィル作用素 $\mathcal{L}_{P}$ と呼ばれ、ハミルトン力学系（エネルギー保存系）で重要になる。
また一方で、クープマン作用素の共役作用素はペロン-フロベニウス作用素 $\mathcal{P}^t$ であり、これを時間微分したものがリウヴィル作用素である。
クープマン固有関数と固有座標
固有値 $\lambda$ に対応する離散クープマン固有関数 $\varphi(\boldsymbol{x})$ は以下を満たす:
\varphi(\boldsymbol{x}_{k+1})
=\mathcal{K}\varphi(\boldsymbol{x}_k)
=\lambda\varphi(\boldsymbol{x}_k)
\tag{21}
連続時間については、固有値を $\mu$ として次の通り:
\frac{d}{dt}\varphi(\boldsymbol{x})
=\mathcal{L}\varphi(\boldsymbol{x})
=\mu\varphi(\boldsymbol{x})
\tag{22}
一般的に、状態空間 $M$ や ダイナミクス（流れ場）$\boldsymbol{F}$ が実数であっても固有値と固有ベクトルは複素数となる。
上の二つの式 $(21)$ と $(22)$ からクープマン固有関数とリー固有関数は固有値は異なるが同じものとなっていることが見て取れる。実際に、定義 $(14)$ を用いて：
\begin{align}
&&\mathcal{K}^t\varphi&=\lambda^t\varphi
\tag{23} \\
&\Rightarrow&
\mathcal{L}\varphi&=\lim_{t\to0}\frac{\mathcal{K}^t\varphi-\varphi}{t}
=\lim_{t\to0}\cfrac{\varphi(\lambda^t-1)}{t}
=\log{(\lambda)}\varphi
\tag{24} \\
&\Rightarrow&
\mu&=\log{(\lambda)} \tag{25} \\
&\Rightarrow&
\lambda^t&=\exp{(\mu t)} \tag{26}
\end{align}
クープマン固有関数と固有ベクトルの関係
離散線形システム $\boldsymbol{x}_{n+1}=\boldsymbol{Ax}_n$ を考える。この時の観測関数の時間発展は $g_{n+1}=\mathcal{K}g_n$ となる。
ここで、行列 $\boldsymbol{A}$ の左固有ベクトル $\boldsymbol{\xi}^\top\boldsymbol{A}=\lambda\boldsymbol{\xi}^\top$ に対応するクープマン固有関数は：
\varphi(\boldsymbol{x})
:=\boldsymbol{\xi}^\top\boldsymbol{x}
\tag{27}
実際に、式 $(21)$ に代入してみると、
\mathcal{K}\varphi(\boldsymbol{x})
=\varphi(\boldsymbol{Ax})
=\boldsymbol{\xi}^\top\boldsymbol{Ax}
=\lambda\boldsymbol{\xi}^\top\boldsymbol{x}
=\lambda\varphi(\boldsymbol{x})
\tag{28}
このように、
- 左固有ベクトルがクープマン固有関数を表し、関数空間 $\mathcal{G}$ における時間不変方向を示す。
- 右固有ベクトル $\boldsymbol{Av}=\lambda\boldsymbol{v}$ は、クープマンモードを表し、状態空間 $M$ における時間不変の方向を与える。
固有値格子
本節では、クープマン固有関数とリー固有関数の格子構造について簡単に記述する。
クープマン作用素の固有値の乗法的性質
クープマン固有関数の集合は新たな固有関数を生成するのに使うことが出来る。2つの固有関数 $\varphi_1(\boldsymbol{x})$ と $\varphi_2(\boldsymbol{x})$ の積 $\varphi_1(\boldsymbol{x})\cdot\varphi_2(\boldsymbol{x})$ もまた固有関数であることがわかる：
\begin{align}
\mathcal{K}
\left(
\varphi_1(\boldsymbol{x})
\varphi_2(\boldsymbol{x})
\right)
&=
\left(
\varphi_1(\boldsymbol{x})
\varphi_2(\boldsymbol{x})
\right)
\circ
\boldsymbol{F}(\boldsymbol{x})
\tag{29} \\
&=
\varphi_1
\left(
\boldsymbol{F}(\boldsymbol{x})
\right)
\varphi_2
\left(
\boldsymbol{F}(\boldsymbol{x})
\right)
\tag{30} \\
&=
\lambda_1\lambda_2
\varphi_1(\boldsymbol{x})
\varphi_2(\boldsymbol{x})
\tag{31}
\end{align}
$\varphi_1$ と $\varphi_2$ の2つの固有値の積で与えられる新しい固有値 $\lambda_1\lambda_2$ に対応する。つまり、クープマン作用素の固有値・固有関数の集合がそれぞれ乗法に閉じていることになる。
リー作用素の固有値の加法性質
式 $(26)$ と $(31)$ の結果を用いれば、
\begin{align}
&&\lambda_1\lambda_2&=e^{\mu_1+\mu_2}
\tag{32} \\
&\Rightarrow&
\mathcal{L}
(
\varphi_1(\boldsymbol{x})
\varphi_2(\boldsymbol{x})
)
&=(\mu_1+\mu_2)
\varphi_1(\boldsymbol{x})
\varphi_2(\boldsymbol{x})
\tag{33}
\end{align}
これはリー作用素の固有値が加法に閉じていることになる。
共役な固有関数から実固有関数の取得
リー作用素の固有関数が複素共役のペア $(\mu,\varphi),\ (\mu ^* ,\varphi ^* )$ として存在する時、実固有関数 $|\varphi|=\sqrt{\varphi\varphi ^*}$ とそれに対応する固有値 $\Re \mu$ が得られる。実固有関数は指数増減のみを表現し、振動成分（=虚部に対応）は持たない。
クープマン固有関数が可換モノイドを成すケース
関数空間 $\mathcal{G}$ が積に閉じている場合、クープマン固有関数空間は可換モノイドを形成する。
可換モノイドについてはこちら
-
乗法に関して可換
$ xy = yx $ -
結合律
$ x \cdot (y \cdot z) = (x \cdot y) \cdot z $ -
単位元の存在
$ x \cdot 1 = 1 \cdot x = x $
故に力学系によっては、全ての固有関数を構成する所謂、基底固有関数集合のようなものが存在する可能性がある。対応するクープマン固有値は先述したように乗法的格子、またはリー固有値の加法的格子を形成する。
観測関数展開
観測関数は固有関数の線型結合として表現できる：
g(\boldsymbol{x})
=\sum_k v_k\varphi_k
\tag{34}\\
ここにクープマン作用素を適用すると：
\begin{align}
\mathcal{K}^tg(\boldsymbol{x})
&=\sum_k \mathcal{K}^tv_k\varphi_k
\tag{35}\\
&=\sum_k v_k\lambda^t_k\varphi_k
\tag{36}
\end{align}
つまり、クープマン作用素が作用（=観測関数の時間変化）すると各固有関数 $\varphi_k$ が各固有値 $\lambda_k^t$ だけスケールすることになる。
クープマン固有関数の等位集合
そもそも等位集合（level set）とは、ある関数が特定の値を取る点の集合のことを指し、例えば関数 $f(\boldsymbol{x})$ に対して、ある定数 $C$ に対する等位集合は：
L_C=\{\boldsymbol{x}\in\mathbb{R}^n\mid f(\boldsymbol{x})=C\}\tag{37}
即ち、等高線のようなものと考えて貰えれば良い。
本節ではクープマン固有関数の等位集合を用いて、システムの構造を明らかにする方法について記述する。特に、固有関数を曲形式で表すことによって、等位集合の性質が直感的に把握し易くなる。
固有関数の極形式表現
固有関数と固有値を極形式で表す：
\begin{align}
\varphi(\boldsymbol{x})&=R(\boldsymbol{x})e^{i\Theta(\boldsymbol{x})}
\tag{38}\\
\lambda&=re^{i\theta}\tag{39}
\end{align}
ここで、$R(\boldsymbol{x})$ は固有関数の大きさを示し、システムの長期的なスケール変化を表す。$\Theta(\boldsymbol{x})$ は固有関数の位相を示し、システムの時間的な進行を表す。
等位集合の性質
続けて、以下2つの集合を考える。
-
等振幅集合
M_\varphi(C) := \{ \boldsymbol{x}\in M: R(\boldsymbol{x})\leq C \} \tag{40}これは $R(\boldsymbol{x})$ が一定以下の集合、つまり等しいエネルギーを持つ領域を表す。
-
等位相集合
A_\varphi(\alpha) := \{ \boldsymbol{x}\in M: \Theta(\boldsymbol{x})\leq\alpha \} \tag{41}これは $\Theta(\boldsymbol{x})$ が一定以下の集合、つまり系の時間進行に対応する等位相曲線となる。
これらの準等位集合は次のような交換法則に従う。$\boldsymbol{x}\in M$ が与えられた時、その次の状態は $\boldsymbol{x}^+=\boldsymbol{F}(\boldsymbol{x})$ と計算でき：
\varphi(\boldsymbol{x}^+)
=
\mathcal{K}\varphi(\boldsymbol{x})
=
\lambda\varphi(\boldsymbol{x})
=
\left(
rR(\boldsymbol{x})
\right)
e^{i(\theta+\Theta(\boldsymbol{x}))}
\tag{42}
つまり、振幅変化はスケール因子 $r$ の乗算、位相変化は一定量 $\theta$ の加算となる。
これをより簡潔に数式を用いると次の通り：
\boldsymbol{F}
(M_\varphi(C))
=
M_\varphi(rC),
\quad
\boldsymbol{F}
(A_\varphi(\alpha))
=
A_\varphi(\alpha+\theta)
\tag{43}
特に、固有値 $\lambda=1$ の場合：
\boldsymbol{F}
(M_\varphi(C))
=
M_\varphi(C),
\quad
\boldsymbol{F}
(A_\varphi(\alpha))
=
A_\varphi(\theta)
\tag{44}
これは、等位集合が時間発展によって変化しない不変集合を形成することを意味する。
クープマン固有関数の推定
クープマン固有関数を推定する方法はいくつか提案されている。ここでは、フーリエ平均によるものと微分方程式によるものを紹介する。
フーリエ平均による推定方法
観測関数の時間発展 $g(\boldsymbol{x}(t))=g(\boldsymbol{F}^t(\boldsymbol{x}))$ を考える。固有値 $e^{i\omega}$ に対応する固有関数は次のようにフーリエ平均を取ることで近似できる：
\tilde{g}_{\omega}(\boldsymbol{x})
:=\lim_{T\to\infty}\frac{1}{T}\int_0^Tg
\left(
\boldsymbol{F}^t(\boldsymbol{x})
\right)
e^{-i\omega t}
\,dt
\tag{45}
これは、$g\left(\boldsymbol{F}^t(\boldsymbol{x})\right)$ の周波数成分 $\omega$ だけを取りだし、時間平均を取る操作になっている。
式 $(42)$ から固有関数は、時間発展を受けても対応する固有値 $e^{i\omega}$ によって単純な変化をする関数であり、周期的なシステムの固有関数の推定に向いている。
偏微分方程式による推定方法
改めて、固有関数を求める永年方程式は：
\mathcal{K}\varphi(\boldsymbol{x})
=
\varphi
\left(
\boldsymbol{F}(\boldsymbol{x})
\right)
=
\lambda\varphi(\boldsymbol{x})
\tag{46}
しかしながらこの合成代数方程式を直接解くのは困難である。そこで、リー作用素に関する永年方程式を考える：
\begin{align}
\mathcal{L}\varphi(\boldsymbol{x})
=
\nabla\varphi(\boldsymbol{x})\cdot
\boldsymbol{f}(\boldsymbol{x},t,\beta)
=
\mu\varphi(\boldsymbol{x})
\tag{47}
\end{align}
この偏微分方程式を解くことで固有関数を近似することが出来る。数値計算的にローラン展開による解法やデータ駆動的に求める解法などがある。
クープマンモード分解と有限表現
ここまで、スカラー値の観測量 $g(\boldsymbol{x})$ に対するクープマン固有関数を考えてきた。本節では複数観測量あるいは（極端に）高次元空間を全て取得する場合を考える。
観測ベクトルの定義
状態 $\boldsymbol{x}\in M$ の観測量を $p$ 個並べてベクトルとして扱う：
\boldsymbol{g}(\boldsymbol{x})
=
\begin{bmatrix}
g_1(\boldsymbol{x}) \\
g_2(\boldsymbol{x}) \\
\vdots \\
g_p(\boldsymbol{x})
\end{bmatrix}
\tag{48}
観測ベクトルの固有関数展開
それぞれの観測量 $g_i(x)$ はクープマン固有関数 $\lbrace\varphi_j(\boldsymbol{x})\rbrace$ を基底として展開できる：
g_i(\boldsymbol{x})
=\sum_{j=1}^\infty v_{ij}\varphi_j(\boldsymbol{x})
\tag{49}
これをベクトル表現に拡張すれば、観測ベクトル $\boldsymbol{g}(\boldsymbol{x})$ は次のように表せられる：
\boldsymbol{g}(\boldsymbol{x})
=
\begin{bmatrix}
g_1(\boldsymbol{x}) \\
g_2(\boldsymbol{x}) \\
\vdots \\
g_p(\boldsymbol{x})
\end{bmatrix}
=\sum_{j=1}^\infty\varphi_j(\boldsymbol{x})\boldsymbol{v}_{j}
\tag{50}
ここで、$\boldsymbol{v}_j$ は固有関数 $\varphi_j(\boldsymbol{x})$ に対応するクープマンモードである。
クープマンモードの計算
ハミルトン力学系では、クープマン作用素が関数空間 $\mathcal{G}(M)\in L_2$ （$L_2$ は自乗可積分集合）上でユニタリ作用素となるため、固有関数が直交規定を成す。
よって、クープマンモード $\boldsymbol{v}_j$ は次のように内積による射影で求めることが出来る：
\boldsymbol{v}_j
=
\begin{bmatrix}
\langle\varphi_j,g_1\rangle \\
\langle\varphi_j,g_2\rangle \\
\vdots \\
\langle\varphi_j,g_p\rangle
\end{bmatrix}
\tag{51}
特に、観測ベクトルが状態そのもの $\boldsymbol{g}(\boldsymbol{x})=\boldsymbol{x}$ である場合、クープマンモードは特定の時間スケールで振動（あるいは成長・減衰）する空間パターンを表す。これらのモード $\boldsymbol{v}$ は動的モード（dynamic modes）としても知られる。
観測ベクトルの時間発展とクープマンモード分解
観測ベクトルの時間発展はクープマン作用素を用いて次のように表せられる：
\boldsymbol{g}(\boldsymbol{x}(t))
=
\mathcal{K}^t\boldsymbol{g}(\boldsymbol{x_0})
\tag{52}
これに式 $(50)$ を適用すると：
\boldsymbol{g}(\boldsymbol{x}_k)
=\sum_{j=1}^\infty
\mathcal{K}^t
\varphi_j(\boldsymbol{x}_0)\boldsymbol{v}_{j}
=\sum_{j=1}^\infty
\lambda^t
\varphi_j(\boldsymbol{x}_0)\boldsymbol{v}_{j}
\tag{53}
この組 $\lbrace(\lambda_j,\ \varphi_j,\ \boldsymbol{v}_j)\rbrace_{j=0}^\infty$ はクープマンモード分解と呼ばれ、後に動的モード分解を介してデータ駆動回帰に繋がる。
クープマンモード分解の近似と動的モード分解
実際にモードを無限個計算するのは不可能なので、支配的なモードを選択することで近似を行う。
クープマンモード分解は理論的手法であり、動的モード分解はデータドリブンな近似手法である。
| 特徴 | クープマンモード分解（KMD） | ダイナミックモード分解（DMD） |
|---|---|---|
| 理論背景 | クープマン作用素に基づく厳密な分解 | クープマン作用素の近似的な固有モード分解 |
| 固有関数 | クープマン固有関数 ( \varphi_j ) | データから推定 |
| 固有値 | クープマン固有値 ( \lambda_j ) | データから推定 |
| モード | クープマンモード ( v_j ) | データから推定 |
| データ | モデルに基づく理論解析 | 時系列データ |
ただし、クープマンモード $\boldsymbol{v}_j$ とクープマン固有関数 $\varphi_j(\boldsymbol{x}_0)$ は異なる数学的対象であり、近似のために異なるアプローチを必要とする。特にクープマン固有関数はより計算困難で、拡張DMDアルゴリズムのような高度な手法が必要となる。
不変固有空間と有限次元モデル
クープマン作用素は本来、無限次元関数空間上で定義されるが、実際には無限次元の進化を捉えるのは困難なので、有限次元の不変部分空間を用いて近似を行う。
ある関数の集合 $\lbrace g_1,g_2,\cdots,g_p\rbrace$ がクープマン不変部分空間を張るとは、任意の線形結合：
g=\sum_j\alpha_jg_j
\tag{54}
が、クープマン作用素 $\mathcal{K}$ を適用した後もこの部分空間に留まることを指す：
\mathcal{K}g
=\sum_j\beta_jg_j
\tag{55}
クープマン作用素を有限個の関数 $\lbrace g_j\rbrace_{j=0}^p$ が張る不変部分空間に制限すると、有限次元の表現行列を得ることが可能となる。この表現行列 $\boldsymbol{K}$ はベクトル空間 $\mathbb{R}^p$ に属し、基底 $g_j(\boldsymbol{x})$ によって求める。（c.f.下図）
有限次元の不変部分空間を構築するには、クープマン固有関数を基底として選択するべきであるが、実際に固有関数を求めるのは困難であるため、近似的な不変部分空間が用いられる。
以下の図は非線形システムに対してクープマン作用素を適用させたイメージ図であるが、$\boldsymbol{y}_k\rightarrow\boldsymbol{x}_k$ の波線は元のダイナミクスの復元を示す。
クープマン埋め込みの例
本節では、固定点を持つ単純な非線形システムをクープマン作用素を用いて解析し、非線形ダイナミクスを線形ダイナミクスとして表現する。具体的には、非線形値を追加してクープマン不変部分空間を構築し、その中で線形表現を求める。
単一の固定点とスローマニホールドを持つ非線形システム
以下のような単一固定点を持つシステムを考える：
\begin{align}
\dot{x}_1&=\mu x_1 \tag{56}\\
\dot{x}_2&=\lambda(x_2-x_1^2)\tag{57}
\end{align}
このシステムは以下の特徴を持つ。
-
$\lambda<\mu<0$ で、$x_2=x_1^2$ に収束するスローマニホールド（ダイナミクスが比較的遅く変化する領域で、挙動が簡単に予測可能な空間）が存在する：
\begin{align} x_1&=Ce^{\mu t}\rightarrow0\quad&\text{if }\mu<0 \tag{58}\\ x_2&\rightarrow x_1^2\quad&\text{if }\lambda<0 \tag{59} \end{align}
さて、このシステムは状態そのものを観測関数 $\boldsymbol{g}(\boldsymbol{x})=(x_1,x_2)^\top$ とするだけでは線形モデルとして扱えない。
そこで、非線形観測量を追加し、$\boldsymbol{g}(\boldsymbol{x})=(x_1,x_2, x_1^2)^\top$ として、3次元のクープマン不変部分空間を構築してみる：
\begin{align}
\begin{bmatrix}
y_1 \\ y_2 \\ y_3
\end{bmatrix}
&=
\begin{bmatrix}
x_1 \\ x_2 \\ x_1^2
\end{bmatrix}
\tag{60} \\
\Rightarrow
\frac{d}{dt}
\begin{bmatrix}
y_1 \\
y_2 \\
y_3
\end{bmatrix}
&=
\frac{d}{dt}
\begin{bmatrix}
x_1 \\ x_2 \\ x_1^2
\end{bmatrix}
\tag{61} \\
&=
\begin{bmatrix}
\mu x_1 \\
\lambda x_2-\lambda x_1^2 \\
2\mu x_1^2
\end{bmatrix}
\tag{62} \\
&=
\begin{bmatrix}
\mu & 0 & 0 \\
0 & \lambda & -\lambda \\
0 & 0 & 2\mu
\end{bmatrix}
\begin{bmatrix}
x_1 \\ x_2 \\ x_1^2
\end{bmatrix}
\tag{63} \\
&=
\begin{bmatrix}
\mu & 0 & 0 \\
0 & \lambda & -\lambda \\
0 & 0 & 2\mu
\end{bmatrix}
\begin{bmatrix}
y_1 \\ y_2 \\ y_3
\end{bmatrix}
\tag{64} \\
\Rightarrow
\cfrac{d}{dt}\,\boldsymbol{y}
&=
\mathcal{K}\boldsymbol{y}
\tag{65}
\end{align}
この3次元空間では、ダイナミクスが完全に線形となり、クープマン作用素の表現行列を記述できる。これを可視化すると、以下のようになる。
可視化デモコード
import numpy as np
import matplotlib.pyplot as plt
# Parameters
lambda_ = -1
mu = -0.05
time = np.linspace(0, 1e4, int(1e6))
# Initial conditions
x1_initial = 2
x2_initial = -1
# Dynamics of x1 and x2
x1 = x1_initial * np.exp(mu * time)
x2 = x1**2 + (x2_initial - x1_initial**2) * np.exp(lambda_ * time)
# Koopman variables y1 = x1, y2 = x2, y3 = x1^2
y1 = x1
y2 = x2
y3 = x1**2
# Plotting
fig = plt.figure(figsize=(8,12))
ax = fig.add_subplot(111, projection='3d')
# Plot the dynamics and the y1-y2 projection
ax.plot(y1, y2, y3, color="black", linewidth=1.0)
ax.plot(y1, y2, -np.ones_like(y3), color="black", linewidth=0.7)
# Plot the y3 = y1^2 surface
Y1, Y2 = np.meshgrid(np.linspace(-2, 2, 100), np.linspace(-1, 4, 100))
ax.plot_surface(Y1, Y2, Y1**2, color="blue", alpha=0.1, rstride=5, cstride=5)
# Plot the y3 = y2 surface
Y1, Y2 = np.meshgrid(np.linspace(-2, 2, 100), np.linspace(0, 4, 100))
ax.plot_surface(Y1, Y2, Y2, color="green", alpha=0.4, rstride=5, cstride=5)
# Plot the y2 = y1^2 surface
Y1, Y3 = np.meshgrid(np.linspace(-2, 2, 100), np.linspace(-1, 4, 100))
ax.plot_surface(Y1, Y1**2, Y3, color="red", alpha=0.1, rstride=5, cstride=5)
# Adding the parabola y = x^2 as a wavy red line in y1-y2 plane]
x_vals = np.linspace(-2, 2, 100)
ax.plot(x_vals, x_vals**2, -np.ones_like(x_vals), color="red", linestyle='--', linewidth=0.7)
# Adding the parabola on the y3 = y2 plane
ax.plot(x_vals, x_vals**2, x_vals**2, color="green", linestyle='--', linewidth=0.7)
# Plot the two points (0, 0, 0) and (0, 0, -1)
ax.scatter(0, 0, 0, color="green", s=30, label="(0, 0, 0)")
ax.scatter(0, 0, -1, color="purple", s=30, label="(0, 0, -1)")
ax.set_xlabel(r'$y_1$')
ax.set_ylabel(r'$y_2$')
ax.set_zlabel(r'$y_3$')
ax.set_xlim(-4,4)
ax.set_ylim(-1,4)
ax.set_zlim(-1,4)
ax.view_init(elev=10, azim=250)
plt.show()
このクープマン埋め込みによって得られた3次元空間には、次の3つの不変多様体が存在する。
- $y_3=y_1^2$（青色）
- 元のシステムにおける不変多様体 $x_1=x_1^2$ に対応
- $y_3=y_2$（緑色）
- 固有値 $\mu, 2\mu$ に対応する固有ベクトルが張るスロー部分空間
- $y_2=y_1^2$（赤色）
- 元のシステムが漸近的に引きつけられる安定多様体
クープマン作用素の固有関数によって定義される固有座標
式 $(27)$ より、クープマン固有関数は左固有ベクトルから求められる。固有値 $\lambda$ に対応する固有ベクトル $\boldsymbol{\xi}$ は次の通り：
\begin{align}
&
\boldsymbol{\xi}^\top\mathcal{K}
=
\lambda\boldsymbol{\xi}^\top
\boldsymbol{O} \tag{66}\\
\Rightarrow
&
\boldsymbol{\xi} =
\begin{bmatrix}
0 \\ 1 \\ -\frac{\lambda}{\lambda-2\mu}
\end{bmatrix}
\tag{67}
\end{align}
これに対応する固有関数 $\varphi_\lambda$ は $\varphi_\lambda=\boldsymbol{\xi}^\top\boldsymbol{x}$ として与えられるので：
\begin{align}
\varphi_\lambda
&=
\begin{bmatrix}
0 && 1 && -\frac{\lambda}{\lambda-2\mu}
\end{bmatrix}
\begin{bmatrix}
x_1 \\ x_2 \\ x_1^2
\end{bmatrix} \\
&=
x_2-\cfrac{\lambda}{\lambda-2\mu}x_1^2,\quad
\end{align}
を得る。
同様にして、固有値 $\mu$ に対しては固有関数 $\varphi_\mu=x_1$ が得られる。
おわりに
次回はこの延長として拡張DMDについて理論紹介を行いたい。