(cache)Koopmanスペクトル解析に関する覚書 - とつとつとしてろうとせず

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

はじめに

動的モード分解（DMD）に関連するデータ解析手法であるKoopmanスペクトル解析（あるいはKoopmanモード分解）の概要についてまとめる。
素人なので理論的な背景など理解していない部分も大いにあるが、なんかすごそう（小並感）な手法であり、今後も発展が期待できる手法と考えている。

元ネタ

とりあえずググってトップに出てくる日本語資料のリンクを貼っておく。今の所、日本語の文献はそんなにない模様（特に流体解析の事例は少ないかも）。
https://www.jstage.jst.go.jp/article/isciesci/61/5/61_175/_pdf
http://www.naoyatakeishi.com/presentation/jsai2017_takeishi_slide.pdf

おすすめは以下のyoutubeの動画。自動字幕ONにすればヒアリング苦手な人（俺）でも何となくわかると思うので是非どうぞ。
www.youtube.com
youtu.be

問題

以降の記述は基本的に上述の文献に基づく。

Koopmanスペクトル解析は一言でいうと「有限次元の非線形問題を無限次元の線形問題に変換して解析する」手法である（なんのこっちゃ）。具体的に見ていこう。

ある状態量 $x$ が時間的に変化する現象に関する問題は次式のような（離散時間）力学系として表現できる。

$\begin{aligned} (1) & x_{k + 1} = f (x_{k}) \end{aligned}$

ここで、 $x_{k}$ は時刻 $t_{k}$ における状態ベクトル、 $f$ は一般に非線形の関数である。したがって、上式の問題を一般的に解く手続きは存在しない。

一方、右辺が線形（かつ時不変）の関数である場合、問題は次式のように表される。

$\begin{aligned} (2) & x_{k + 1} = A x_{k} \end{aligned}$

ここで、 $A$ は適当な行列である。DMDの記事で見たように、この差分方程式は一般的に解くことができる。そのため、解けない非線形の方程式を近似的にでも解ける線形の方程式に変換する、というのが非線形問題に対する基本的なアプローチになる。

Koopman作用素

さて、次式を満たす作用素 $K$ をKoopman作用素（あるいは合成作用素）と呼ぶ。

$\begin{aligned} (3) & K g (x_{k}) = g (f (x_{k})) = (g \circ f) (x_{k}) \end{aligned}$

ここで、 $g$ は状態量に対する観測量（observable）と呼ばれる任意の関数であり、これ自体は非線形であってもよい。定義から分かる通り、Koopman作用素は観測量と考えている力学系の右辺の関数の合成関数を作るものである。このKoopman作用素により、状態量 $x$ に関する問題を観測量 $g$ に関する問題として書き換えることができる。

$\begin{aligned} (4) & g (x_{k + 1}) = K g (x_{k}) \end{aligned}$

つまり、Koopman作用素は観測量の下での時間発展を表す作用素であることがわかる。 $K$ は関数空間に対して働くので一般にその次元は無限であるが、しかし作用素としては線形である。実際、 $a, b$ をスカラーとして次式が成り立つ。

$\begin{aligned} (5) & K (a g_{1} + b g_{2}) (x) = a K g_{1} (x) + b K g_{2} (x) \end{aligned}$

逆にいえば、あえて次元を無限次元にまで上げることによって問題を線形なものに変換できたともいえる。というわけで、Koopman作用素のアイディアは「元の状態量による表現では非線形であった問題を、適当な観測量をうまく選んでやることで線形の問題に変換する」ものといえる。

先に離散時間系について考えたが、連続時間系についても同様に考えることができる。

$\begin{aligned} (6) & \frac{d}{d t} x = f (x) \end{aligned}$

この場合、初期値 $x_{0}$ からスタートして $f$ によって描かれる状態空間上の軌道 $F^{t} (x_{0})$ を考える（つまり、上式を積分した式）。

$\begin{aligned} (7) & F^{t} (x_{0}) = x (t) = x_{0} + \int_{0}^{t} f (x (τ)) d τ \end{aligned}$

このとき、連続時間系の場合のKoopman作用素は次式で定義される。

$\begin{aligned} (8) & K^{t} g (x_{0}) = g (F^{t} (x_{0})) \end{aligned}$

線形作用素の固有値と固有関数

線形作用素については行列と同じように固有値と固有ベクトル（固有関数）を考えることができる。俺自身関数解析に明るくないのでイメージ重視の理解だが、この概念は微分方程式を解く際に有用であり、Sturm–Liouville理論などが有名な例だと思う。

たとえば、単純な例として次式を考える。

$\begin{aligned} (9) & \frac{d}{d t} x (t) = λ x (t) \end{aligned}$

このとき、微分作用素 $\frac{d}{d t}$ の固有値は任意の実数 $λ \in R$ であり、固有関数はこの微分方程式の解 $e^{λ t}$ である。一般に作用素の固有値は離散的なものと連続的なものを両方持ち得るそうだが、ここでは基本的に離散的な固有値の話を取り上げる。

さて、Koopman作用素も線形作用素であるため、固有値と固有関数を考えることができる。

$\begin{aligned} (10) & K ϕ_{j} (x) = λ_{j} ϕ_{j} (x) \end{aligned}$

ただし、実際の計算ではKoopman作用素を陽な形で書き下すことは（少数の例外を除いて）通常できない。そこで、代わりに上式を満たすKoopman固有関数を探すことが課題となり、研究が進められているようだ（候補関数の組み合わせで近似する方法やカーネル法を使う方法、ニューラルネットワークなど機械学習を援用する方法など）。

簡単な例として、有限次元の線形力学系を考える。

$\begin{aligned} (11) & x_{k + 1} = A x_{k} \end{aligned}$

このとき、この系におけるKoopman固有関数は $ϕ_{j} (x) = ⟨ x, w_{j} ⟩$ となる。ただし、 $⟨, ⟩$ は適当な内積、 $w_{j}$ は $A$ の随伴固有ベクトル（実数なら左固有ベクトル）である（ $A^{*} w_{j} = λ_{j}^{*} w_{j}$ ）。
なぜなら、随伴作用素（随伴行列）の定義より、 $⟨ A x, y ⟩ = ⟨ x, A^{*} y ⟩$ なので、

$\begin{aligned} (12) & \begin{aligned} K ϕ_{j} (x_{k}) & = ϕ_{j} (A x_{k}) = ⟨ A x_{k}, w_{j} ⟩ = ⟨ x_{k}, A^{*} w_{j} ⟩ \\ = ⟨ x_{k}, λ_{j}^{*} w_{j} ⟩ = λ_{j}^{*} ⟨ x_{k}, w_{j} ⟩ \\ = λ_{j}^{*} ϕ_{j} (x_{k}) \end{aligned} \end{aligned}$

となるからだ。

また、Koopman固有関数は次のような非常に素晴らしい性質を持つ。新しい状態量をKoopman固有関数を用いて $z_{j} = ϕ_{j} (x)$ とすると、

$\begin{aligned} (13) & \begin{aligned} z_{j} [k + 1] = ϕ_{j} (x_{k + 1}) = ϕ_{j} (f (x_{k})) = (K ϕ_{j}) (x_{k}) = λ_{j} ϕ_{j} (x_{k}) = λ_{j} z_{j} [k] \end{aligned} \end{aligned}$

となるので、まとめてベクトル形で書けば次式を得る。

$\begin{aligned} (14) & \begin{aligned} z_{k + 1} = Λ z_{k} \end{aligned} \end{aligned}$

すなわち、Koopman固有関数による座標変換によって非線形力学系を線形化することができる（これはTaylor級数による局所的な線形化ではなく、大域的なものである点に注意）。この意味でもKoopman固有関数は重要であるといえる。なお、座標変換による大域的な線形化は制御工学の分野で研究されているようだ。

ちなみに、連続時間系であれば、Koopman固有関数を求めるための線形な偏微分方程式を導くことができる（参考文献はこれとかこれ）。連続時間系の場合、

$\begin{aligned} (15) & \begin{aligned} \frac{d}{d t} ϕ (x) = K ϕ (x) = λ ϕ (x) \end{aligned} \end{aligned}$

であるので、微分の連鎖律より、

$\begin{aligned} (16) & \begin{aligned} \frac{d}{d t} ϕ (x) = \nabla_{x} ϕ (x) \cdot \dot{x} = \nabla_{x} ϕ \cdot f (x) = λ ϕ (x) \\ ∴ \nabla_{x} ϕ \cdot f (x) = λ ϕ (x) \end{aligned} \end{aligned}$

という偏微分方程式を得る。制御系など連続力学系の陽な表現が分かっている場合は使えそうだ（離散的な時系列データから回帰法によって求めることもできるらしい）。

Koopman固有関数が簡単には見つからないことの例を挙げる（上述の動画から引用）。次の1次元力学系を考える。

$\begin{aligned} (17) & \dot{x} = f (x) = x^{2} \end{aligned}$

上式は非線形といっても変数分離で積分できてしまうので問題にならない、というのには目をつぶってもらうことにして、上式を観測量をうまく選ぶことで線形化できないかと考えてみる。素朴に $g (x) = x^{2}$ としてみよう。状態空間表現に書き直すと、

$\begin{aligned} (18) & \frac{d}{d t} [\begin{matrix} g_{1} \\ g_{2} \end{matrix}] = [\begin{matrix} 0 & 1 \\ ? & ? \end{matrix}] [\begin{matrix} x \\ x^{2} \end{matrix}] \end{aligned}$

という形で書くことができれば線形化できるということだ。しかし、微分の連鎖律より、

$\begin{aligned} (19) & {\dot{g}}_{2} = \frac{d g_{2}}{d x} \cdot \frac{d x}{d t} = 2 x \cdot \dot{x} = 2 x \cdot x^{2} = 2 x^{3} \end{aligned}$

となるので、 $g_{1}, g_{2}$ だけでは閉じた式にならないことがわかる。じゃあ、新しく $g_{3} (x) = x^{3}$ とおいてみたらどうか。しかし、全く同様に ${\dot{g}}_{3} = 3 x^{4}$ となるので、この連鎖は終わらない。この手順を無限回繰り返せば次のような無限次元の線形系を得る。

$\begin{aligned} (20) & \frac{d}{d t} [\begin{matrix} g_{1} \\ g_{2} \\ g_{3} \\ ⋮ \end{matrix}] = [\begin{matrix} 0 & 1 & 0 & 0 & \dots \\ 0 & 0 & 2 & 0 & \dots \\ 0 & 0 & 0 & 3 & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}] [\begin{matrix} x \\ x^{2} \\ x^{3} \\ ⋮ \end{matrix}] \end{aligned}$

上式はKoopman作用素が次元を無限次元に上げることで問題を線形化する好例といえるだろう。しかし、実際の計算上、無限次元は扱えないわけなので、上式の連鎖はどこかで区切りよく終わってもらうか、せめて高次項ほど影響が小さくなって途中で打ち切っても誤差を小さくできるような形が望ましい。というわけで、最適なKoopman固有関数を探すことが課題となる。

ちなみに、上式は先述の偏微分方程式から固有関数を求めることができる。1次元の場合、

$\begin{aligned} (21) & \frac{d ϕ}{d x} f (x) = λ ϕ (x) \end{aligned}$

となるので、

$\begin{aligned} (22) & \frac{d ϕ}{d x} = λ \frac{ϕ (x)}{f (x)} \end{aligned}$

という常微分方程式を得る。これを積分すると、

$\begin{aligned} (23) & ϕ (x) = ϕ_{0} \exp (\int_{x_{0}}^{x} \frac{λ}{f (s)} d s) = ϕ_{0} e^{- λ / x} \end{aligned}$

したがって、定数倍の違いを除いて求める固有関数は $ϕ (x) = e^{- 1 / x}$ であることがわかる。実際、時間微分を計算してみると、

$\begin{aligned} (24) & \dot{ϕ} = x^{- 2} e^{- 1 / x} \cdot \dot{x} = ϕ (x) \end{aligned}$

となり、観測量を適切に、すなわち、固有関数に選ぶことで非線形問題を線形化できることがわかる。もちろん、上記の偏微分方程式を用いる方法は万能ではなく、数値的なデータでしか表現できない離散時間系に対しては固有関数を得ることは一般に難しい課題である。

Koopmanモード分解（KMD）

Koopman固有関数に基づくモード分解手法がKMDである。いま、観測量 $g$ がKoopman固有関数によって次のように展開可能であるとする。

$\begin{aligned} (25) & \begin{aligned} g_{i} (x_{k}) = \sum_{j = 1}^{\infty} ϕ_{j} (x_{k}) v_{i j} \end{aligned} \end{aligned}$

ここで、 $v_{i j}$ は展開係数である。つまり、観測量をKoopman固有関数によって張られる部分空間の中から選ぶということだ。このとき、

$\begin{aligned} (26) & \begin{aligned} ϕ_{j} (x_{k}) = ϕ_{j} (f (x_{k - 1})) = (K ϕ_{j}) (x_{k - 1}) = λ_{j} ϕ_{j} (x_{k - 1}) = λ_{j}^{k} ϕ_{j} (x_{0}) \end{aligned} \end{aligned}$

より、観測量のベクトルをまとめて次のように書くことができる。

$\begin{aligned} (27) & \begin{aligned} g (x_{k}) = \sum_{j = 1}^{\infty} λ_{j}^{k} ϕ_{j} (x_{0}) v_{j} \end{aligned} \end{aligned}$

$v_{j}$ は展開係数をまとめたベクトルであり、Koopmanモードと呼ばれる。実際には $ϕ_{j} (x_{0}) v_{j}$ の形で求まる。

DMDとの関係

Koopmanモード分解を実際に求めるアルゴリズムは、実はDMDそのものである。というのも、DMDは観測量として $g (x) = x$ 、すなわち、状態量そのものを使ったKMDと解釈できるからだ。また、データとなる観測量としてKoopman固有関数を取ることができれば、DMDによって求まる時間発展行列はKoopman作用素の有限次元部分空間における近似となる。というわけで、DMDはKMDのサブジャンルに位置する手法といえる。

まとめ

Koopmanスペクトル解析に関する概要についてまとめた。面倒なんで具体的な計算例だったり実行例だったりは上述の動画などを見てくだせえ（ｵｲ）。実際、本手法のエッセンスについては動画の方がわかりやすいので。