定在波のDMD - とつとつとしてろうとせず

注意

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

概要

復習がてらDMDの資料作ってたらハマったのでメモ。以前調べたときに理解してたつもりだったのだが、すっかり忘れていたという。やっぱり自分で手を動かしてみないとダメですね。

この記事でいいたいこと：Tu(2014)の3.2.2節をよーく読んで理解しておくように。

DMDモードと周波数

動的モード分解（DMD）は複雑な時系列データから特徴的な構造（DMDモード）を取り出すデータ解析手法の一つなのであった。同様の手法に固有直交分解（POD）があるが、DMDの長所はPODと違い、空間モードと時間周波数が同時に求まることである。
というわけで、PODとDMDを比較する例としてよく紹介されている計算をやってみよう。ちなみに実行環境はJupyterLabである（notebookよりこっちのが好きかも。RISEが使えればなおよし）。
次のような関数のモードを求めることを考える。
$\begin{align} \begin{aligned} f(x, t) &= x e^{j t} + \cos(x) e^{2 j t} + \tanh(x) e^{3 j t} \\ & -\pi \leq x \leq \pi \\ & 0 \leq t \leq 6\pi \end{aligned} \end{align}$ ここで、 $j = \sqrt{-1}$ である（Pythonの記法に準拠）。入力データがあからさまにモード分解されているわけだが、問題はこれらのモードをデータだけから復元できるか、ということである。
計算のPythonコードも載っけようかと思ったけど整理が面倒だったので省略。勉強したい人はこことかここを見てほしい（実行結果付きでわかりやすい）。
上記の関数を3Dプロットすると下図のようになる（実部だけ）。
f:id:iqujack-lequina:20181202111321p:plain
これにDMDを実行し、空間モードを取り出すと下図のようになる（これも実部だけ。データとしては複素数である点に注意）。横に並べたのはPODモード（という名の左特異ベクトル）である。時間平均は引いてないけど、$\int_{0}^{2n\pi} e^{j m \pi t} dt = 0$（$ n, m $は自然数）なので問題ないはず。なお、特異値分解の打ち切り次数は3とした。
f:id:iqujack-lequina:20181202115446p:plain
DMDはきちんと空間モードが抽出できているのに対し、PODモードは元の関数形を再現できていないのがわかる。
連続時間に変換した固有値を見てみると以下のようになる。

[ -8.29622288e-15+1.j  -4.90814069e-15+3.j  -8.42442059e-15+2.j]

きちんと虚部が設定した数値になっている（実部も数値誤差程度でほぼ0だ）。一方、PODは常に実数の特異値しか求まらないため、当然周波数の情報はない。PODとDMDの違いがよくわかるシンプルな例である。
一応、時間項の実部のプロットもつけておく。山の数が周波数に一致している。
f:id:iqujack-lequina:20181202175548p:plain

実数のデータ

さて、ここで次のような疑問を抱く人もいるかもしれない（ぶっちゃけ俺のことだけど）。

「サンプルの関数は$e^{j \omega t}$じゃなくて単に$\sin(\omega t)$とか$\cos(\omega t)$でよくね？　実際のデータって大体実数なわけだし」

確かに、お説ご尤も。では実際にやってみよう。
上のデータから実部だけ取り出したデータに対して同様にDMDをやってみる。
空間モードを見てみると、さっきと同様の結果が得られた（符号の正負は係数で調整できるのでどちらでもよい）。
f:id:iqujack-lequina:20181202175304p:plain
しかし、固有値を見てみるとどうだろう、見事に実部だけになってしまった。

[-0.04743168+0.j -0.43210939+0.j -0.19058688+0.j]

固有値の虚部が周波数を表すのだから、これじゃ振動モードが取り出せていないじゃん、となる。DMD使えねーな！
だがちょっと待ってほしい。ちゃんとTu(2014)にこの問題についての説明がある（個人的にはこっちのバージョンの説明の方がわかりやすい）。Standing wave、つまり定在波の問題だ。
DMDではモードの周波数は時間発展行列の複素固有値の虚部として現れる。逆にいえば、固有値が実数のみの場合、周波数の情報が出てこないので、DMDは時間的な変動をうまく捉えることができない。最も単純な例は $x_{k} = \sin(\omega t_{k})$ という1次元の時系列データに対してDMDを適用する場合である（上の例でいえば、ある特定の空間点の時間変動のみを見ている場合に相当）。このとき、データ行列は1次元の行ベクトルになるが、そうすると、時間発展行列は $1 \times 1$ の行列、つまり、ただのスカラーになってしまう。当然、固有値も1つのスカラーとなり、元データが実数なので固有値も実数である。先程の例も同様で、各行が定在波になっているため、固有値がすべて実数となっている。
では、こうしたデータにはDMDは適用できないのかというと、そうでもない。元のデータが1次元であっても、1つの周波数につきモードの空間は2次元である必要があるが、逆にいえば、不足している次元を補ってやればよいということになる。
これは次のような簡単な手続きで行うことができる。すなわち、元データから次式のようなデータ行列を作成すればよい。
$\begin{align} \begin{aligned} \mathbf{Z} = \begin{bmatrix} \mathbf{X}_{0:m-1} \\ \mathbf{X}_{1:m} \end{bmatrix} \end{aligned} \end{align}$ ここで、$\mathbf{X}_{k:l}$は時刻$t_{k}$から時刻$t_{l}$までのデータを列方向に並べた行列である。つまり、DMDのアルゴリズムで使う2個の行列を行方向に並べた行列ということだ。Pythonで書くとこんな感じ。

# shift-stacked data
Z = np.vstack((X[:,:-1],X[:,1:]))

これにより、線形独立な行ベクトルが2倍に増えるので、共役複素数によって周波数を表現できるようになる。あとはこの新しいデータ行列$\mathbf{Z}$を改めて2つに分け、DMDを実行すればよい。ただし、結果の行数が2倍になっているので、最終的には半分の行までのデータを見るようにする。
結果は下図の通りである。今度は打ち切り次数を倍の6に設定しているため、同じモードが2つずつあるが、やはり空間モードは正しく計算できる*1。
f:id:iqujack-lequina:20181202175457p:plain
固有値を見ると、今度はきちんと共役複素数のペアが得られた（見やすさ優先で虚部だけ表示しているが、実部はほぼ0である）。

[ 3. -3.  1. -1.  2. -2.]

時間遅れ座標

上記の「時間をずらしたデータを新しいデータと見なす」手法は、最初に知ったときは1次元データを無理矢理DMDの定式化に落とし込むためにデータを水増ししているようにしか思えなかった。しかし、このような考え方は「時間遅れ座標」として知られており、カオスな力学系のアトラクタを再構成する場合などに用いられているらしい（参考文献）。これは2階以上の高階の線形常微分方程式を解くために、微分階数に応じた変数を追加して高次元化（ベクトル化）し、1階の常微分方程式として解くというおなじみの手順にも似ているかもしれない（時刻をずらしたデータを使うことが微分に対応している）。また、Koopman固有関数を求める場合にも使用されている（参考文献）。そういった意味で、単なるデータの水増しではなく、積極的な意味のある手法のようである。
時間遅れ座標を用いた解析にも興味はあるので、やる気と時間があったら別の機会にまとめたい（この論文も読んだけど、よくわからないんだよなぁ）。