動的モード分解に関する覚書 - とつとつとしてろうとせず

注意

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

はじめに

流体解析などの分野でデータ解析や縮約モデリングの手法として用いられる動的モード分解の概要についてまとめる。近年、PODと並んで注目されている？手法であり、主に海外で盛んに研究されているようだ（国内での流行り具合はどうなんだろ）。実際には分野や状況に応じてPODと使い分けていくことになるのかな。
関連手法であるKoopmanモード分解については別記事でまとめる予定。

元ネタ

DMDの概要と実行例が載っている素晴らしいサイト。
基本的にはこちらを見てもらえばいいと思うが、PODの場合と違ってネットに日本語の解説があまりないようなので、勉強として自分でもまとめておく。

POD

まずはPOD（固有直交分解）についておさらいすると、時系列データを並べて作ったデータ行列に対して主成分分析を行い、データから特徴的な構造（固有モード）を抽出するのがPODなのだった（名前の由来はPODモードが正規直交基底となるため）。
具体的な計算としてはデータ行列を特異値分解することで行われる。特異値の大きなモードを支配的なモードと解釈し、主要なモードのみを抜き出すことで元のデータからデータ容量を削減した縮約モデル（ROM）を作り、低コストなシミュレーションを行うこともできる。
また、時系列データだけでなく、パラメータを変えて得られたデータに対しても適用できるので、最初に代表的なパラメータに対してデータをサンプリングしておけば、パラメータが動く範囲内に対してはROMによって高速なシミュレーションができる（データを補間するイメージ）。

DMD

PODによる流れ場のモード分解では、時間と空間の相関を同時に計算するため、時空間のモードがうまく分離されず、特に時間的な振る舞いが複数のモードに分散して抽出されてしまう場合がある（系を特徴付ける1つの周波数を再現するために複数の自明でないモードの組み合わせが必要になる）。
動的モード分解（Dynamic Mode Decomposition, DMD）では、時系列データの時間発展に注目してモード分解を行うことで、固有モードとそれに対応する周波数を得ることができるので、あるモードが時間的に成長するのか減衰するのか、または周期的に振動するのかといった、いわゆる安定性に関する情報を得ることもできる。

以下、具体的なアルゴリズムについてまとめる。定式化としてはKrylov部分空間法に基づくものと特異値分解に基づくものがあるが、数学的な解析には前者が使われ、数値計算には後者が用いられるようなので、ここでは特異値分解に基づくアルゴリズムについて述べる。

離散力学系の時間発展

PODと同様、データとしては各時刻におけるスナップショットを1次元のベクトル（状態ベクトル）としてまとめたものを用いる。データ数を $m$ 個とすると、このデータベクトル$\mathbf{x}_k$の時間発展は次式のように表される。

$\begin{align} \mathbf{x}_{k+1} = \mathbf{f} \paren{ \mathbf{x}_{k} } \quad \paren{ k = 0, 1, 2, \cdots , m-1 } \end{align}$

ここで、$\mathbf{f}$は系の時間発展を表す関数のベクトルであり、一般に非線形で場合によっては未知である。上式は数学的には離散力学系と呼ばれる定式化になっており、ある時刻の状態が次の時刻にどのような状態に移るのかを表す（連続力学系の場合は上式は時間に関する一階の微分方程式になる）。
DMDはこの時間発展をデータから推定するという発想に基づいており、データ駆動型（data-driven ）で方程式によらない（equation free）手法といわれる。つまり、物理モデルが存在しない未知の現象やそもそもモデル化が困難な経済学的な現象（株価の変動とか？）など、$\mathbf{f}$が具体的に分からない現象に対しても、データだけを使ってその振る舞いを分析・予測しようというアプローチを取っている。

さて、上式をまとめて行列で書くと、次式のようになる。

$\begin{align} \mathbf{X}' = \mathbf{F} \paren{ \mathbf{X} } \end{align}$

ここで、データ行列 $\mathbf{X} \in \mathbb{R}^{n \times m-1}$ と $\mathbf{X}' \in \mathbb{R}^{n \times m-1}$ は次式で定義される。

$\begin{align} \mathbf{X} = \begin{bmatrix} \vert & \vert & & \vert \\ \mathbf{x}_{0} & \mathbf{x}_{1} & \cdots & \mathbf{x}_{m-2} \\ \vert & \vert & & \vert \end{bmatrix} \\ \mathbf{X}' = \begin{bmatrix} \vert & \vert & & \vert \\ \mathbf{x}_{1} & \mathbf{x}_{2} & \cdots & \mathbf{x}_{m-1} \\ \vert & \vert & & \vert \end{bmatrix} \end{align}$

つまり、上式は元の $m$ 個のデータを並べた行列から、時刻を1つずつずらした行列になっており、$\mathbf{X}$と$\mathbf{X}'$の相関関係から時間発展を求めていくことになる。そこで、DMDでは上式を次式で近似する。

$\begin{align} \mathbf{X}' \approx \mathbf{A} \mathbf{X} \end{align}$

すなわち、未知の非線形関数を線形の行列（線形写像）で近似する。元々線形ならともかく、非線形の場合にもこういうことをして上手くいくのかという疑問は当然であって、実際上手くいかないこともあるようだ。ただし、ある状況下では非線形の時間発展を線形近似しても上手くいくことが知られており、この辺りの話はKoopman作用素によるモード分解と関連している（詳細については別記事で触れる）。

力学系の安定性

ちょっと脱線して力学系の話についてまとめる。門外漢のため、俺も詳しくはないのだが、この辺の概念がDMDのアルゴリズムや結果を理解するのに必要になるっぽいので、最低限イメージだけでも掴めるようにメモっておく。

連続力学系

まずは連続力学系から考える。一般の状態ベクトルの時間発展は次の微分方程式で表される。

$\begin{align} \dif{}{t} \mathbf{x} = \mathbf{f} \paren{ \mathbf{x}} \end{align}$

右辺が$t$に陽に依存する場合は非自律系（非自励系）、そうでない場合は自律系（自励系）という。以降、簡単のため自律系を考える。
上述のように $\mathbf{f}$ は一般には非線形であるため、そのまま上式を解くことはできない。そこで、速度が0になる点（平衡点）の周りで上式を線形化し、近似的にその振る舞いを調べるということが行われる。したがって、適当な線形化の下で次式を得る。

$\begin{align} \dif{}{t} \mathbf{x} = \mathbf{A} \mathbf{x} \end{align}$

上式は線形の微分方程式であり、容易に次のような解を得ることができる。

$\begin{align} \mathbf{x}(t) = e^{\mathbf{A} t} \mathbf{x}(0) \end{align}$

ここで、 $e^{\mathbf{A} t}$ は行列指数関数と呼ばれ、具体的には $\mathbf{A}$ に関する固有分解を $\mathbf{A} = \mathbf{U}\boldsymbol{\Lambda}\mathbf{U}^{-1}$ として、次式で与えられる。

$\begin{align} \begin{aligned} e^{\mathbf{A} t} &= \mathbf{U} e^{\boldsymbol{\Lambda} t} \mathbf{U}^{-1} \\ e^{\boldsymbol{\Lambda} t} &= \operatorname{diag} \paren{ e^{\lambda_{j} t} } \end{aligned} \end{align}$

つまり、 $\mathbf{A}$ が対角化可能であれば、行列指数関数は $\mathbf{A}$ の固有値と固有値ベクトルから計算できる（対角化ができなくてもJordan標準形を用いることで同様に計算可能）。

指数関数の性質とEulerの公式から、 $\lambda_{j} = \alpha_{j} + i \omega_{j}$ とすれば、

$\begin{align} e^{\lambda_{j} t} = e^{\alpha_{j} t} \paren{ \cos\paren{\omega_{j} t} + i \sin\paren{\omega_{j} t} } \end{align}$

であり、すべての固有値の実部が負（ $\alpha_{j} < 0$ ）なら系は安定であり、そうでなければ不安定である。したがって、線形系（平衡点）の安定性は行列の固有値から判定できる。

離散力学系

離散力学系に対しても同様に考えることができる。この場合、平衡点に対応するのは $\mathbf{f} \paren{ \mathbf{x} } = \mathbf{x}$ を満たす固定点（あるいは不動点）である。これは差分商（というかEuler陽解法？）

$\begin{align} \dif{}{t} \mathbf{x} \approx \frac{ \mathbf{x}_{k+1} - x_{k} }{ \Delta t } = \frac{ \mathbf{f} \paren{ \mathbf{x}_{k} } - x_{k} }{ \Delta t } = 0 \end{align}$

を考えると分かりやすいかもしれない。

さて、適当な固定点で線形化すると次式を得る。

$\begin{align} \mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k} \end{align}$

上式は差分方程式であり、時刻 $k$ の解 $\mathbf{x}_{k}$ は初期値から繰り返し行列を作用させることで次のように得られる。

$\begin{align} \begin{aligned} \mathbf{x}_{1} &= \mathbf{A} \mathbf{x}_{0} \\ \mathbf{x}_{2} &= \mathbf{A} \mathbf{x}_{1} = \mathbf{A}^{2} \mathbf{x}_{0} \\ \mathbf{x}_{3} &= \mathbf{A} \mathbf{x}_{2} = \mathbf{A}^{3} \mathbf{x}_{0} \\ & \vdots \\ \therefore \mathbf{x}_{k} &= \mathbf{A}^{k} \mathbf{x}_{0} \end{aligned} \end{align}$

行列のべき乗は対角化によって簡単に計算することができ、

$\begin{align} \mathbf{A}^{k} = \mathbf{U} \mathbf{\Lambda}^{k} \mathbf{U}^{-1} \end{align}$

となるので、離散力学系の場合、固定点はすべての固有値の絶対値が1より小さければ安定であり（ $\abs{\lambda_{j}} < 1$ ）、そうでなければ不安定である。

連続系と離散系の関係

連続系と離散系の対応関係は以下のようになる。区別のため、それぞれの時間発展行列を $\mathbf{A}_{c}, \ \mathbf{A}_{d}$ とすると、離散時間が $t_{k} = k \Delta t$ で与えられることを思い出せば、

$\begin{align} \begin{aligned} \mathbf{x} \paren{(k+1) \Delta t} &= e^{\mathbf{A}_{c} \paren{(k+1)\Delta t}} \mathbf{x} \paren{0} \\ &= e^{\mathbf{A}_{c} \Delta t } e^{\mathbf{A}_{c} \paren{k \Delta t}} \mathbf{x} \paren{0} \\ &= e^{\mathbf{A}_{c} \Delta t } \mathbf{x} \paren{k \Delta t} \end{aligned} \end{align}$

であるから、

$\begin{align} \mathbf{A}_{d} = e^{\mathbf{A}_{c} \Delta t} \end{align}$

となることが分かる。また、連続系の対角化行列を $\mathbf{U}_{c}$ とすると、

$\begin{align} \begin{aligned} \mathbf{A}_{d} \mathbf{U}_{c} &= e^{\mathbf{A}_{c} \Delta t} \mathbf{U}_{c} \\ &= \paren{ \mathbf{U}_{c} e^{\boldsymbol{\Lambda}_{c} \Delta t} \mathbf{U}_{c}^{-1} } \mathbf{U}_{c} \\ &= \mathbf{U}_{c} e^{\boldsymbol{\Lambda}_{c} \Delta t} \end{aligned} \end{align}$

となることから、連続系と離散系の固有ベクトルは等しく（ $\mathbf{U}_{d} = \mathbf{U}_{c}$ ）、離散系の固有値は $\boldsymbol{\Lambda}_{d} = e^{\boldsymbol{\Lambda}_{c} \Delta t}$ で与えられることが分かる（逆にいえば、離散系の固有値から連続系の固有値を $\lambda_{c} = \frac{\ln \lambda_{d}}{\Delta t}$ で計算できる）。

また、 $\abs{e^{\lambda \Delta t}} = \abs{e^{\paren{\alpha + i \omega}\Delta t} } = \abs{e^{\alpha \Delta t}} \abs{e^{i \omega \Delta t}}$ より、連続系が安定（ $\alpha < 0$ ）なら対応する離散系も安定であることがわかる（ $\Delta t > 0 ,\ \abs{e^{i \omega \Delta t}} = 1$ なので $\abs{e^{\alpha \Delta t} } < 1$ となる）。

時間発展行列の計算

ともあれ、次なる問題はこの行列$\mathbf{A}$をどう求めるかだ。最初に考えるのは$\mathbf{X}$の逆行列を使って$\mathbf{A} = \mathbf{X}' \mathbf{X}^{-1}$とすればいいじゃんというものだが、そもそもの作り方からしてデータ行列は基本的に非正方行列であるし、かりに正方行列であったとしても正則性（逆行列の存在）は保証されない。そこで、近似的に上式が成り立つような行列を次式の条件から求めるという方針になる（数式の書き方これであってんのかな……？）。

$\begin{align} \mathbf{A} \triangleq \operatorname*{argmin}_{\mathbf{A}} \| \mathbf{X}' - \mathbf{A} \mathbf{X} \|_{F} \end{align}$

ここで、$\| \cdot \|_F$はFrobeniusノルム（行列ノルム）である。上式を満たす行列はMoore-Penroseの擬似逆行列$\mathbf{X}^{\dagger}$を使って次式のように得られる。

$\begin{align} \mathbf{A} \approx \mathbf{X}' \mathbf{X}^{\dagger} \end{align}$

$\mathbf{X}$の特異値分解が求まっていれば、擬似逆行列は次式で計算できる。

$\begin{align} &\mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{*} \\ \Longrightarrow \quad & \mathbf{A} = \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \mathbf{U}_{r}^{*} \end{align}$

ただし、$^*$は複素共役転置である（DMDの計算では固有値や固有ベクトルが複素数になる）。また、下付きの $r$ は特異値分解の時点でランク$r$の低ランク近似を行うことを意味する（特異値の大きな順に $r$ 個の要素まで残してそれ以降は打ち切る）。

後は行列$\mathbf{A}$の固有値・固有ベクトルを求めてやれば時間発展の周波数と固有モードを得ることができる。

ただし、$\mathbf{A} \in \mathbb{C}^{n \times n}$より、行列$\mathbf{A}$のサイズが非常に大きいため、そのまま固有値問題を解くことは得策ではない。そのため、実際の計算では次式のようにPOD基底 $\mathbf{U}_{r}$ に射影した行列$\tilde{\mathbf{A}}$を求め、こちらに対して固有値問題を解くことで計算コストを小さくする（相似変換であるため、固有値は同じになる）。

$\begin{align} &\tilde{\mathbf{A}} = \mathbf{U}_{r}^{*} \mathbf{A} \mathbf{U}_{r} = \mathbf{U}_{r}^{*} \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \\ \Longrightarrow \quad & \tilde{\mathbf{A}} \mathbf{W} = \mathbf{W} \boldsymbol{\Lambda} \end{align}$

したがって、元の行列 $\mathbf{A}$ の固有モード（DMDモード）が次式のように得られる。

$\begin{align} \begin{aligned} \mathbf{A} \boldsymbol{\Phi} &= \boldsymbol{\Phi} \boldsymbol{\Lambda} \\ \boldsymbol{\Phi} &= \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \mathbf{W} \end{aligned} \end{align}$

DMDモードを上式で求める手法をexact DMDと呼ぶ。一方、元々の相似変換の関係からDMDモードを $\boldsymbol{\Phi} = \mathbf{U}_{r} \mathbf{W}$ から求める手法はprojected DMDと呼ぶ（DMDが提案されたオリジナルの論文ではこちらを使っている）。exact DMDは演算に $\mathbf{U}_{r} = \mathbf{X} \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1}$ ではなく、 $\mathbf{U}'_{r} = \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1}$ を使っている点に注意。
Tu(2014)によれば、projected DMDは $\mathbf{X}$ のデータしか使っておらず、もはや元の行列 $\mathbf{A}$ の固有ベクトルにはなっていないが、exact DMDは $\mathbf{X}$ と $\mathbf{X}'$ の両方のデータを使っており、 $\mathbf{A}$ の固有ベクトルになっている。そのため、exact DMDの方が厳密であり、自然な手法ということらしい（実際、今はexact DMDの方が主流っぽい）。

DMDモードを用いた予測

さて、固有モードが求まれば解の再構成が可能になる。離散力学系 $\mathbf{x}_{k+1} = \mathbf{A} \mathbf{x}_{k}$ の解は固有分解 $\mathbf{A} = \boldsymbol{\Phi} \boldsymbol{\Lambda} \boldsymbol{\Phi}^{\dagger}$ を用いて次式で与えられるのだった（ $\boldsymbol{\Phi}$ が正則とは限らないので擬似逆行列になる）。

$\begin{align} \mathbf{x}_{k} = \boldsymbol{\Phi} e^{\boldsymbol{\Lambda}^{k}} \boldsymbol{\Phi}^{\dagger} \mathbf{x}_{0} \end{align}$

離散系と連続系の対応関係から、 $\boldsymbol{\Omega} = \ln \boldsymbol{\Lambda} \left/ \Delta t \right.$ とおけば、上式を任意の時間に対する解として表現することができる（離散的なデータから連続な時間の解が求まる！）。

$\begin{align} \mathbf{x}(t) = \boldsymbol{\Phi} e^{\boldsymbol{\Omega t}} \mathbf{b} \end{align}$

ただし、 $\mathbf{b} = \boldsymbol{\Phi}^{\dagger} \mathbf{x}_{0}$ とおいた。これは次のようなDMDモードによるモード展開を表している点に注意する（ $\boldsymbol{\xi}$ は適当な空間座標）。

$\begin{align} \mathbf{x}(t) = \sum_{j=0}^{r} b_{j}(t) \varphi_{j}\paren{\boldsymbol{\xi}} = \sum_{j=0}^{r} b_{j}(0) e^{\omega_{j} t} \varphi_{j}\paren{\boldsymbol{\xi}} \end{align}$

PODと異なり、こちらは時間発展を計算するだけなら展開係数の計算は要らないのが特徴だ。原理的には周期的な現象のデータを1周期分だけ取っておけば、それ以降の振る舞いを何周期分でも計算で予測することができる（とはいえ、実際には長時間予測には課題が多いようだ）。
また、DMDモードは直交基底ではないということもPODと異なる点だ。これは一見悪いことのように思えるが、直交性の条件に縛られないことによって柔軟な表現力があるともいえる（主成分分析と独立成分分析の関係に似ているかも）。ROMシミュレーションを考える場合は双対基底（左固有ベクトル）を使うことで展開係数も簡易に計算できる（Wang&Wei(2017)）。

計算手順

DMDのアルゴリズムをまとめると、以下のようになる。

時系列データからデータ行列$\mathbf{X}$と$\mathbf{X}'$を作成
$\mathbf{X}$の特異値分解 $\mathbf{X} = \mathbf{U}_{r} \boldsymbol{\Sigma}_{r} \mathbf{V}_{r}^{*}$ を計算し、低ランク近似を行う
POD基底に射影した時間発展行列を $\tilde{\mathbf{A}} = \mathbf{U}_{r}^{*} \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1}$ で求める
$\tilde{\mathbf{A}}$ の固有値問題 $\tilde{\mathbf{A}} \mathbf{W} = \mathbf{W} \boldsymbol{\Lambda}$ を解く
元の行列の固有モード（DMDモード）を $\boldsymbol{\Phi} = \mathbf{X}' \mathbf{V}_{r} \boldsymbol{\Sigma}_{r}^{-1} \mathbf{W}$ で計算
任意の時刻に対する解は $\mathbf{x}(t) = \boldsymbol{\Phi} e^{\boldsymbol{\Omega t}} \boldsymbol{\Phi}^{\dagger} \mathbf{x}_{0}$ で計算できる（ $\boldsymbol{\Omega} = \ln \boldsymbol{\Lambda} \left/ \Delta t \right.$ ）

まとめ

DMDの概要についてまとめた。本当ならpythonでテスト計算もやるべきなんだけど、力尽きたので別の記事でやることにする（実行例は上述のブログを見てください）。
個人的には離散系と連続系の対応関係について確認できてよかった（論文とか読んでて固有値のlogを取ってるのを見て訳がわからなかったので）。
やっぱり力学系とか現代制御とかをもっと勉強しないとかなぁ。