動的モード分解(Dynamic Mode Decomposition)を実装してみた

はじめに

こんにちは．かわつあんです．最近，私の研究の一つで動的モード分解(Dynamic Mode Decomposition, 以下DMD)というものを使っています．DMDは元来，流体解析などの高次元数をもつシミュレーションデータの次元削減を目的として発展してきた，データ駆動の解析の一種です．最近現れた解析で，理論的な研究がすすみ，私の流体解析分野や非線形力学系，そして私の研究分野のプラズマ物理におけるシミュレーション・実験への応用が徐々に行われてきています．

DMDの利点は教師となるような関数をひつようとしないデータ駆動的な計算で，主要モードとその周波数および減衰率まで計算できる点です．(フーリエ解析はsin波を，ウェーブレット解析ではマザーウェーブレットを教師としていますよね)．また，フーリエ解析と異なり時間分解能と周波数の最小幅が対応していません．一方で，DMDでは変調された構造や周波数が時間変化するような系ではその情報を取り出せないというデメリットもあります．

今回の記事では，DMDとそのアルゴリズムの説明とDMDの実装をいます．なお，コーディングは以下のサイトを参考にしています． www.pyrunner.com

私の理解に誤りがありましたら，是非コメントでご指摘等してくださると幸いです．

動的モード分解(DMD)

まずは普通のDMDを実装します． DMDの歴史や思想などは極力省略して，アルゴリズムの説明に重きを置いていますことをご留意ください．

各時間 $t_k, ( k=0,1, \cdots, m )$ 毎に得られた長さ $n$ の列ベクトルを $\bf x_k$ とします．これを普通に並べると

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

となると思います. DMDでは，はじめに次のように並べた２つ行列 $\bf X, Y$ を考えます．

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

これはすなわち， $\bf X$ は最後の時間のデータ以外を並べた行列であり， $\bf Y$ は最初の時間のデータ以外を並べた行列になっており， $\bf X$ と $\bf Y$ は時間的に１ステップだけズレています．実験やシミュレーションで得られるデータは非線形的で動的なものがほとんどですが，DMDでは動的なデータについて線形であると近似します（詳くは [J. H. Tu (2014)]などを参考）．

すなわち

$\mathbf{Y = AX}$

と書けると考えます． $A$ はある行列であり線形作用素です． $A$ を求めるためには

$\mathrm{min} \bf || AX - Y||_F$

となる $A$ を求めればよいはず．ここで， $\bf ||P||_F$ はフロベニウスノルム( $\bf P$ はなにかしらの行列)です．一般に $\bf X, Y$ は正方行列ではないので，一般的な行列の逆行列として定義される擬似逆行列を考えます．擬似逆行列 $\bf X^+$ とすると

$\bf A \approx YX^+$

として作用素 $\bf A$ を定められます．DMDでは作用素 $\bf A$ に対する固有値方程式を解き，得られた固有ベクトルからDMDモードを計算し，固有値からDMDモードの周波数および減衰率を求めます．

一般的なDMDでは擬似逆行列 $\bf X^+$ を求めるために特異値分解を行います．特異値分解は以下のように書けます．

$\bf X = U \Sigma V^H$

ここで， $\bf X$ は $n \times m$ 行列でしたが， $\bf U$ は $n \times n$ のユニタリー行列， $\bf \Sigma$ は $n \times m$ の対角行列， $\bf V$ は $m \times m$ のユニタリー行列です．また， $\bf V^{H}$ は $\bf V$ のエルミート共役です．エルミート共役とはある行列を転置して複素共役をとったもので，実数行列における転置行列に対応しています．ユニタリー行列はある行列とエルミート共役の積が単位行列になる場合，その行列はユニタリー行列であるといい（つまり $\bf U^{H} U = U U^{H} = E$ ， $\bf E$ は単位行列），実数行列の直行行列に対応しています．

以上のように特異値分解があるすると擬似逆行列 $\bf X^+$ は

$\bf X^+ = V \Sigma^{-1} U^H$

とかけるので， $A$ の射影行列 $\tilde{A}$ を

$\bf \tilde{A} = U^H A U = U^H YX^+U = U^HYV\Sigma^{-1}$

として求め， $\bf \tilde{A}$ についての固有値方程式

$\bf \tilde{A} w = \Lambda w$

を解きます．得られた $\bf \tilde{A}$ の固有値は $\bf A$ の固有値でもあるため，DMD アルゴリズムの目的である $\bf A$ の固有値方程式が解けたことになります．

standard-DMDの場合，固有ベクトルは以下の通り計算すると，

$\begin{align} \mathbf{\tilde{A}w} &=\mathbf{w\Lambda} \\ \mathbf{U^HAUw} &=\mathbf{w\Lambda} \\ \mathbf{A(Uw)} &=\mathbf{(Uw)\Lambda} \\ \end{align}$

となるため，固有ベクトルは $\bf Uw$ となります．

一方で，現在一般的に使われているexact-DMDの場合， $A$ に対応する固有値ベクトルは，以下のように書き直されます．

$\bf \Phi = YV\Sigma^{-1}w$

確認のため， $\bf A \Phi$ を計算すると，

$\begin{align} \mathbf{A \Phi} &=\mathbf{A Y V\Sigma^{-1}w } \\ &=\mathbf{(YV\Sigma^{-1}U^H )Y V\Sigma^{-1}w} \\ &=\mathbf{YV\Sigma^{-1}\tilde{A} w}\\ &=\mathbf{YV\Sigma^{-1}w\Lambda}\\ &=\mathbf{\Phi\Lambda}\\ \end{align}$

となり，確かに $\bf \Phi$ が $\bf A$ の固有ベクトルであることができました．以上がDMDの基本的なアルゴリズムですが，基本的なアルゴリズムにもいろいろあるそうなので，参考まで．

DMDを実装

では，先ほど説明したDMDを実装してみましょう．まずテストデータを作成します．

import numpy as np
import scipy.linalg                          #scipyの線形数学モジュール．　略記は非推奨らしい
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid

#時空間を定義
x = np.linspace(-10,10,60)
t = np.linspace(0,20,80)
Xm, Tm = np.meshgrid(x, t)

#テストデータ生成
D1 = 5 * (1/np.cosh(Xm/2)) * np.tanh(Xm/2) * np.exp((0.8j)*Tm) #主要モード
D2 = 0.2 * np.sin(2 * Xm) * np.exp(2j * Tm)                    #サブモード
D3 = 0.1 * np.random.randn(*Xm.shape)                          #ノイズ
D = (D1 + D2 + D3)
fig = plt.figure(figsize = (12,3))

#データをプロット
#Primary mode
sp1 = plt.subplot(1,4,1)
sp1.contourf(x, t, D1, cmap="jet", levels = np.linspace(-3,3,50))
sp1.set_aspect(1)
sp1.set_xlabel("x")
sp1.set_ylabel("t")
sp1.set_title("primary mode")
sp1.set_yticks([0,5,10,15,20])
#Secondary mode
sp2 = plt.subplot(1,4,2)
sp2.contourf(x, t, D2, 50, cmap="jet", levels = np.linspace(-3,3,50))
sp2.set_aspect(1)
sp2.set_xlabel("x")
sp2.set_ylabel("t")
sp2.set_title("secondary mode")
sp2.set_yticks([0,5,10,15,20])
#noise
sp3 = plt.subplot(1,4,3)
sp3.contourf(x, t, D3, 50, cmap="jet", levels = np.linspace(-3,3,50))
sp3.set_aspect(1)
sp3.set_xlabel("x")
sp3.set_ylabel("t")
sp3.set_title("noise")
sp3.set_yticks([0,5,10,15,20])
#created data
sp4 = plt.subplot(1,4,4)
cont=sp4.contourf(x, t, D, 50, cmap="jet", levels = np.linspace(-3,3,50))
sp4.set_aspect(1)
sp4.set_xlabel("x")
sp4.set_ylabel("t")
sp4.set_title("created data")
sp4.set_yticks([0,5,10,15,20])
plt.tight_layout()

実行結果

f:id:kawats1:20200120011449p:plain — テストデータ

このテストデータに対してDMDモードを計算してみます．先ほどのプログラムに追記の形で以下のように実装します．

def dmd(data, truncation = None):
    #step1　データから行列を生成
    X = data[: ,1:]
    Y = data[:,:-1]

    #step2　特異値分解
    U, Sigma, V_H = scipy.linalg.svd(X, False) #scipy.linalg.svdでは，特異値は特異値のリストとして出力される
    
    #step3　射影行列計算　（ランクの切り捨ても行う．）特異値のランクの切り捨ては引数truncationで指定
    r = len(Sigma) if truncation is None else truncation
    
    U = U[:,:r]
    V = V_H.conj().T[:,:r]
    U_H = U.conj().T
    Sigma = np.diag(Sigma)[:r,:r] #特異値のリストを対角行列に．
    Sigma_inv = scipy.linalg.inv(Sigma)
    #step4 A_tildeの固有値問題を解く．
    A_tilde = np.dot(np.dot(np.dot(U_H, Y),V), Sigma_inv)
    Lambda, W = scipy.linalg.eig(A_tilde) 
    #step5 固有値に対応するAに固有ベクトルを求める
    Phi = np.dot(np.dot(np.dot(Y,V), Sigma_inv), W )
    return Lambda, Phi

# 
# 空間パターン
eigen_val, modes = dmd(D,truncation = 3)
fig1 = plt.figure()
plt.title("Spatial Pattern")
plt.plot(x, modes)
plt.legend(["mode1","mode2","mode3"])
plt.xlabel("x")
fig1.savefig("SP")

# 時空間パターン
omega=np.log(eigen_val)/0.1              # 複素角振動数
fig2 = plt.figure(figsize = (6,3))
plt.suptitle("Spatialtemporal Pattern")
sp5=plt.subplot(1,2,1)
sp5.contourf(x, t, np.outer(W_[:,0],np.exp(omega[0] * t)).T, 50, cmap="jet" ) 
sp5.set_aspect(1)
sp5.set_xlabel("x")
sp5.set_ylabel("t")
sp5.set_yticks([0,5,10,15,20])
sp6=plt.subplot(1,2,2)
sp6.contourf(x, t, np.outer(W_[:,1]*(W_[:,1]),np.exp(omega[1] * t)).T,50, cmap="jet" ) 
sp6.set_aspect(1)
sp6.set_xlabel("x")
sp6.set_ylabel("t")
sp6.set_yticks([0,5,10,15,20])
plt.tight_layout()

実行結果

f:id:kawats1:20200120005602p:plain — DMDモードの空間パターン

f:id:kawats1:20200120005635p:plain — DMDモードの時空間パターン

無事，DMDを使って２つのモードの抽出に成功しました．（ノイズについてはあまり関係ないので，今回は省略しています．）なお，固有値方程式を解いて得られた固有値 $\lambda_j$ から，複素角振動数 $\omega_j$ を次の通り計算しています．

$\omega_j = \mathrm{log}(\lambda_j)/ \Delta t$

もちろん実周波数 $f_j$ ，成長率 $\nu_j$ はそれぞれ以下のように計算できます．

$f_j = \mathfrak{Im}(\omega_j)/2\pi,\quad \nu_j = \mathfrak{Re}(\omega_j)/2\pi$

ここで注意が一つあります．実は今回実装したのはDMD解析の中でも空間パターンの基底関数やその時間発展を求めただけです．この他にもDMDモードの振幅を求めたり，DMDモードの重ね合わせで元のモードを再現できるかどうかのvalidationなども必要ですが，今回は割愛しています．また，詳しい理論などにも触れておらず申し訳ありません．

そこらへんの内容は需要があればまとめますが，基本的には原著論文などを参考にしてください

＝＞参考

以上のコードはGistでまとめていますので参考にどうぞ．(すこしMarkdownが文字化け？してますが．．．)

gist58c04e83be338c5786317c114fdde62e