初めに

表題の

POD（Proper Orthogonal Decomposition, Lumley1967, Berkooz1993）
DMD（Dynamic Mode Decomposition, Schmid2010, Tu2013, Tu2014 (arxiv), Schmid2022）

は似ているように感じるのだが，，，恥を忍んで正直に白状すると，よく分かっていない．（幾つかのvariantsがあることは知っているのだが，取り敢えず標準の）DMDのことを知っておきたい．

参考に，例えば，平2011a，平2011b．周辺に，左海+2014，Wu+2019．個人的に興味のある話題に，Brunton+2022．まあ，この辺りは，Brunton, Kutz, Taira, Herrmann, McKeonの庭だろうか．

Methods

Data

方法論を理解する，という枠組みでは，具体的な例を持ち出すのは頭が固くなる要因になり得るのだが，，，しかし具体的なイメージがあった方が理解し易いというのも事実である．取り敢えず，過去に作ったコードを使う．

$\text{Re}$	$\text{Results}$
$ 20$
$ 100$
$ 200$
$2,000$

よく見るように$\text{Re} = 100$を対象に，$10 \sim 15 \ [\mathrm{s}]$のperiodic flowを考えてみる．

ところで，どうでも良いことだが，離散化パラメータ$h = 1 \times 10^{-2} \ [\text{m}]$，$\tau = 1 \times 10^{-3} \ [\text{s}]$程度で刻んだこの結果は，$\text{Re} = 100$程度では，Norberg1993, Sharma&Eswaran2004, Cheng+2007, Sen+2011, Mashhadi+2021などと比較して，悪くない結果であった．AMRを用いたLBMの数値解（Fakhari&Lee2014）とも定性的には合致している．$\text{Re} = 1,000$程度でも，円柱周りの流れ（e.g., Behr+1995, Durante+2017）に合致する様子がある．本当はもっとblockageを抑えたりするのだろうが，既にやり尽くされた問題であろう．

POD

== 独り言ここから ==
私の理解が正しいかは置いておいて，いまの私の理解は，この動画と合致している．この動画が正解という訳でもないし，この動画を基礎に勉強した訳でもない．ただ，後からこの動画を見つけて，ああ，確かに僕もそういう理解なんだよね，と感じた．だから，この動画でも言及されていない要素までを私が理解し切れていない可能性は十分ある．
== 独り言ここまで ==

各時刻での渦度$\boldsymbol{w}$を縦ベクトルに並び替え，その時系列を行方向に並べた行列$\boldsymbol{W}$を用意する．

\begin{align}
    \mathbb{K}^{M \times N}
    \ni
    \boldsymbol{W}
    &= \left[ 
        \boldsymbol{w}_0, \boldsymbol{w}_1, \ldots, \boldsymbol{w}_{N-1}
    \right]
\end{align}

ここで，$M$は空間の自由度，$N$はタイムステップ数である．えいやっという数字だが，空間を$0.01 \ [\mathrm{m}]$で刻んだ$M \sim 400 \times 100 \sim 40 \ \mathrm{K}$程度の節点で，$0.01 \ [\mathrm{s}]$ごとのスナップショットを5秒間だけ用意した（$N \sim 500$）．ところで，平均値を除いておき，ゼロ周りまで平行移動させておく．これに対して，SVDを行う．

\begin{align}
    \boldsymbol{W}
    &= \boldsymbol{U}
        \boldsymbol{S}
        \boldsymbol{V}^{*}
     = \sum_{r}
        s_{r}
        \boldsymbol{u}_{r}
        \boldsymbol{v}^{*}_{r}
\end{align}

ここで，$\boldsymbol{U}, \ \boldsymbol{V}$はユニタリ，$\boldsymbol{V}^{*}$は随伴である．$\boldsymbol{S}$は特異値を成分に持つ対角行列である．得られた特異値と，その累積量の関係を見ると，上位の幾つかのモードの内に，90％以上の情報が閉じ込められていることが分かる．

一つひとつのモードは，次のような見た目をしている．界隈では，空間的特徴量のことをモードと呼び，時間的特徴量のことをダイナミクスと呼ぶようである．個人的には，空間モード・時間モードと呼んだ方がしっくりくるのだけれど．．．ところで，色々なものに対してSVDを行ったことがあるが，往々にしてペアを見る，何故なのだろう？

0番目は平均場である．事前に各スナップショットの平均値を除いていても，このような平均場が残るのは予想に反していたが，考えれば当然のようにも感じる．水平流速についても同じことをやれば色々分かる，というメモ（ゼロ周りか否か）．

また，パラメータの変動による解の変化を特徴付ける空間を探すような使い方もできる．昔，大学院の講義で学んだが，興味深い内容であった．なお，当該の講義では，これを次数低減化法とか低減化基底法とか呼んでいたが，最近の論文を幾つか読むと，そういうのはPOD-Galerkin projection methodとか呼ばれている印象がある．文化だろうね．

どうでも良いことだが，mode 5, 6辺りを見ると，いつもムキムキの腹筋を思い浮かべてしまう．

DMD

いま，渦度（に着目しているから渦度の話をするが何でも良い）自体はNS式を解くことで観測してきた．しかし，この渦度の時間発展を支配する方程式が存在するはずである（NS式に回転を作用させたり？．．．取り敢えず未知とする）．即ち，

\begin{align}
    \frac{d}{d t} \boldsymbol{w}
    &= f \left(
        \boldsymbol{w}
    \right)
\end{align}

或いは，離散的に書けば，

\begin{align}
    \frac{\boldsymbol{w}_{n+1} - \boldsymbol{w}_{n}}{\tau}
    &= f \left(
        \boldsymbol{w}_{n}
    \right)
    \\
    \therefore
    \boldsymbol{w}_{n+1}
    &= \boldsymbol{w}_{n}
        + \tau f \left(
        \boldsymbol{w}_{n}
    \right)
    \\
    &= F \left(
        \boldsymbol{w}_{n}
    \right)
\end{align}

DMDでは，この$F \left( \cdot \right)$を観測データから推定する．この枠組みのことを，Brunton2018は，equation-freeと，インパクトのある言葉で宣伝している．

具体的な手続きとして，PODのときと同様に，各時刻での物理量を縦ベクトルに並び替え，その時系列データを行方向に並べた行列を用意する．但し，時間にズレを作り，2つの行列を作る．

\begin{align}
    \mathbb{K}^{M \times N}
    \ni
    \boldsymbol{W}_{0}
    &:= \boldsymbol{W}_{0:N-1}
     = \left[ 
        \boldsymbol{w}_0, \boldsymbol{w}_1, \ldots, \boldsymbol{w}_{N-1}
    \right]
    \\
    \mathbb{K}^{M \times N}
    \ni
    \boldsymbol{W}_{1}
    &:= \boldsymbol{W}_{1:N}
     = \left[ 
        \boldsymbol{w}_1, \boldsymbol{w}_2, \ldots, \boldsymbol{w}_{N}
    \right]
\end{align}

これらを用いて，先に登場した離散化した時間発展を書くと，

\begin{align}
    \boldsymbol{W}_{1}
    &= F \left(
        \boldsymbol{W}_{0}
    \right)
\end{align}

しかし，手がかり無しに考えるのは難しい．えいやっと線形写像で時間発展が記述されるのだと考える（線形で良いのか甚だ疑問であるが，そういう方針を取り，工夫として，Koopman作用素を噛ませるとかあるらしい）．

\begin{align}
    \boldsymbol{W}_{1}
    &= F \left(
        \boldsymbol{W}_{0}
    \right)
    \\
    &\simeq \boldsymbol{A}
        \boldsymbol{W}_{0}
        \quad \left(
            \boldsymbol{A}   \in \mathbb{K}^{M \times M}, \
            \boldsymbol{W}_0, \boldsymbol{W}_1 \in \mathbb{K}^{M \times N}
        \right)
\end{align}

即ち，$F \left( \cdot \right)$を求める問題から，$\boldsymbol{A}$を求める問題へと置き換えた．$\boldsymbol{A}$を求めるには，

\begin{align}
    \boldsymbol{W}_{1}
    &= \boldsymbol{A} \boldsymbol{W}_{0}
    \\
    &= \boldsymbol{A} \boldsymbol{U} \boldsymbol{S} \boldsymbol{V}^{*}
    \\
    \therefore
    \boldsymbol{A}
    &= \boldsymbol{W}_{1} \boldsymbol{W}_{0}^{+}
    \\
    &= \boldsymbol{W}_{1}
        \left(
            \boldsymbol{U}_{0} \boldsymbol{S}_{0} \boldsymbol{V}_{0}^{*}
        \right)^{+}
    \\
    &= \boldsymbol{W}_{1}
        \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1} \boldsymbol{U}_{0}^{*}
\end{align}

である（$\boldsymbol{M}^{+}$は pseudoinverse of the Moore-Penrose type, メモ書き）．但し，$\boldsymbol{A} \in \mathbb{K}^{M \times M}$はクソデカなので（$M$は空間の自由度），具体的な計算はやりたくない．一方，線形近似しているので，$\boldsymbol{A}$が全てであることも事実だ．したがって，$\boldsymbol{A}$を求める代わりに，$\boldsymbol{A}$の固有値や固有ベクトルを求める問題へ置き換える（したがって，DMDとは，$\boldsymbol{A}$を直接計算すること無しに$\boldsymbol{A}$を求めることである，と解釈している）．

この上で，$\boldsymbol{A}$を$\boldsymbol{U}_{0}$に写す（と，Tu2013やBrunton2018は言っているが，どういうことか？．．．良く分かっていない）．

\begin{align}
    \mathbb{K}^{N \times N}
    \ni
    \tilde{\boldsymbol{A}}
    :&= \boldsymbol{U}_{0}^{*} \boldsymbol{A} \boldsymbol{U}_{0}
    \\
    &= \boldsymbol{U}_{0}^{*} \boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1}
\end{align}

暫くの間，分からないと思っていたのだが，次のようなことかもしれないと思い始めた．元の（線形近似した）系は，

\begin{align}
    \boldsymbol{W}_{1}
    &= \boldsymbol{A}
        \boldsymbol{W}_{0}
\end{align}

である．ここで，$\boldsymbol{W}$たちを変換する（というか，これを吐き出すような潜在的な空間を考える，と言った方が良いのだろうか）．

\begin{align}
    \boldsymbol{W}_{0}
    &= \boldsymbol{U}_{0} \tilde{\boldsymbol{W}}_{0}
    ,\
    \boldsymbol{W}_{1}
     = \boldsymbol{U}_{0} \tilde{\boldsymbol{W}}_{1}
        \quad \left(
            \tilde{\boldsymbol{W}}_{0}, \ \tilde{\boldsymbol{W}}_{1} \in \mathbb{K}^{N \times N}
        \right)
    \\
\end{align}

すると，

\begin{align}
    \boldsymbol{W}_{1}
    &= \boldsymbol{A}
        \boldsymbol{W}_{0}
    \\
    \boldsymbol{U}_{0} \tilde{\boldsymbol{W}}_{1}
    &= \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{W}}_{0}
    \\
    \tilde{\boldsymbol{W}}_{1}
    &= \boldsymbol{U}_{0}^{*} \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{W}}_{0}
    \\
    &= \tilde{\boldsymbol{A}} \tilde{\boldsymbol{W}}_{0}
    \quad \left(
        \tilde{\boldsymbol{A}} := \boldsymbol{U}_{0}^{*} \boldsymbol{A} \boldsymbol{U}_{0}
    \right)
\end{align}

であり，低次元のダイナミクスが存在するように見える．Koopmanのような雰囲気を感じる．一応，筋は通る気がするので，こういう解釈で進む．

なお，$\tilde{\boldsymbol{A}}$の定義の書き換えには，

\begin{align}
    \boldsymbol{W}_{1}
    &= \boldsymbol{A} \boldsymbol{W}_{0}
     = \boldsymbol{A} \boldsymbol{U}_{0} \boldsymbol{S}_{0} \boldsymbol{V}_{0}^{*}
    \\
    \Leftrightarrow
    \boldsymbol{A} \boldsymbol{U}
    &= \boldsymbol{W}_{1} (\boldsymbol{S}_{0} \boldsymbol{V}_{0}^{*})^{+}
     = \boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1}
\end{align}

を使った．$\boldsymbol{A} \in \mathbb{K}^{M \times M}$よりも$\tilde{\boldsymbol{A}} \in \mathbb{K}^{N \times N}$の方が小さく，嬉しい．$\tilde{\boldsymbol{A}} \in \mathbb{K}^{N \times N}$の固有値分解を行う．

\begin{align}
    \tilde{\boldsymbol{A}}
    &= \tilde{\boldsymbol{L}}
        \tilde{\boldsymbol{\Lambda}}
        \tilde{\boldsymbol{L}}^{-1}
    \Leftrightarrow
    \tilde{\boldsymbol{A}}
        \tilde{\boldsymbol{L}}
     = \tilde{\boldsymbol{L}}
        \tilde{\boldsymbol{\Lambda}}
\end{align}

これは，$\boldsymbol{A}$のスペクトル$\lbrace \boldsymbol{\Lambda}, \boldsymbol{L} \rbrace$の代わりに，$\tilde{\boldsymbol{A}}$のスペクトル$\lbrace \tilde{\boldsymbol{\Lambda}}, \tilde{\boldsymbol{L}} \rbrace$を調べよう，というマインドだと思うのだが，両者の固有値や固有ベクトルは一致しているのだろうか？という疑問が残っている（一致することが示せるらしいのだが，，，次元も違うのに一体どうなっているのか？cf. こちら, こちら）．ところで，$\tilde{\boldsymbol{A}}$は非対称であろうから，$\tilde{\boldsymbol{\Lambda}}$は複素数だろう．

先の固有値分解の式に，左から$\boldsymbol{A} \boldsymbol{U}_{0}$をぶつける．

\begin{align}
    \tilde{\boldsymbol{A}} \tilde{\boldsymbol{L}}
    &= \tilde{\boldsymbol{L}} \tilde{\boldsymbol{\Lambda}}
    \\
    \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{A}} \tilde{\boldsymbol{L}}
    &= \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{L}} \tilde{\boldsymbol{\Lambda}}
\end{align}

分けて計算して，

\begin{align}
    \text{(LHS)}
    &= \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{A}} \tilde{\boldsymbol{L}}
    \\
    &= \boldsymbol{A} \boldsymbol{U}_{0} (\boldsymbol{U}_{0}^{*} \boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1}) \tilde{\boldsymbol{L}}
    \\
    &= \boldsymbol{A} \boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1} \tilde{\boldsymbol{L}}
    \\
    \text{(RHS)}
    &= \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{L}} \tilde{\boldsymbol{\Lambda}}
    \\
    &= (\boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1}) \tilde{\boldsymbol{L}} \tilde{\boldsymbol{\Lambda}}
    \\
\end{align}

より，

\begin{align}
    \tilde{\boldsymbol{A}} \tilde{\boldsymbol{L}}
    &= \tilde{\boldsymbol{L}} \tilde{\boldsymbol{\Lambda}}
    \\
    \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{A}} \tilde{\boldsymbol{L}}
    &= \boldsymbol{A} \boldsymbol{U}_{0} \tilde{\boldsymbol{L}} \tilde{\boldsymbol{\Lambda}}
    \\
    \boldsymbol{A} (\boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1} \tilde{\boldsymbol{L}})
    &= (\boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1} \tilde{\boldsymbol{L}}) \tilde{\boldsymbol{\Lambda}}
\end{align}

これを，$\boldsymbol{A}$の固有値分解と解釈する．即ち，

\begin{align}
    \boldsymbol{L}
    := \boldsymbol{W}_{1} \boldsymbol{V}_{0} \boldsymbol{S}_{0}^{-1} \tilde{\boldsymbol{L}}
\end{align}

が固有ベクトルである，と云う．これをDMDモードと呼ぶ．これは exact DMD と呼ばれているらしい．一方，$\boldsymbol{L} := \boldsymbol{U} \tilde{\boldsymbol{L}}$と書くこともできて，これは projected DMD と呼ばれているらしい．

結局，DMDがやっているのは，

$F (\cdot)$を探す問題を，$\boldsymbol{A}$を探す問題に置き換え，
$\boldsymbol{A}$を探す問題を，$\tilde{\boldsymbol{A}}$を探す問題に置き換え，
$\tilde{\boldsymbol{A}}$を探す問題を，$\lbrace \tilde{\boldsymbol{\Lambda}}, \tilde{\boldsymbol{L}} \rbrace$を探す問題に置き換え，
$\tilde{\boldsymbol{\Lambda}}$を$\boldsymbol{\Lambda}$として採用することで$\tilde{\boldsymbol{L}}$から$\boldsymbol{L}$を回復する

ものであると解釈した．誤っているかもしれない．

Results

PODについては一定の経験がある．しかし，DMDについてきちんと勉強するのは初めてだったので，数式の議論と数値の議論は分けておきたかった．

Re = 100 (periodic)

PODのときと同じように，$\mathrm{Re} = 100$のperiodic flowに対してこれを適用する．

初めに計算するのは$\boldsymbol{W}_{0}$に対するSVDであるが，これはPODと何ら変わらないので省略する．

次に登場する具体的な数値計算は，$\tilde{\boldsymbol{A}}$のスペクトルである．上から10個だけ見てみる．なお，前述の通り，$\tilde{\boldsymbol{A}}$は非対称なので，スペクトルは複素数である．

$\lambda_6 = 1 + 0 i$がstationaryであり，空間場を見てもPODの第0次モードと同様の平均場としての表現がある．

$\lambda_9 = 0.9 + 0 i$は単位円より内側にあり，エネルギーの減衰を示唆する．実際，$\exp (\log (\lambda_9) / \Delta t \cdot t)$が減衰している．しかし，僅か1秒以内に減衰し切っており，空間場も極めて高周波の，ノイズのような成分であろうと考えられる．何故なら，十分に発達したあとのperiodicな流れを見ており，殆ど減衰していない．

Re = 100 (developing regime)

冒頭の動画からも分かるように，流れがperiodicな状態で安定するまでに数秒かかる．今回の計算では，初期に与えたノイズが抜け，渦が発達し始める$3 \sim 8 \ [\mathrm{s}]$の状態に対してPOD，DMDを適用してみる．

POD

ネットで探しても中々見ない結果を得られておーという感じ．mode 1-2は少しずつ振幅を増しており，mode 7-9辺りはノイズが消えてゆく過程のように見える．

DMD

少々見辛いが，，，減衰するモードと増幅するモードが重ね合わせられている．mode 4と5が全く逆さまのように見える．

Re = 2,000 (developing regime)

先ほどと同じく，渦が発達し始める早い時間（$5 \sim 10 \ [\mathrm{s}]$）の間にPOD，DMDを適用する．

POD

0番目の平均場が振動している．都合の良い定数が見当たらないくらいゆらゆらしているということだろう．また，再度mode 7-9周辺には細かい波が乗っている．しかし，これはもうノイズではなく，本当に細かい波たちなのではないか，と思っている

DMD

多くが減衰してゆき，僅かな成分だけがstable，或いは増幅している．少し期待外れだったが，動画を見返すと確かに落ち着いてゆく時間帯かもしれない．

Re = 2,000 (noise-free)

これまでは，渦の発達を促進するために，初期条件に擾乱を与えていた．静かな流れが渦を作るまで待ってみて，場の様相が激しく変化する状態に対してPOD，DMDを適用してみる．

0~5 [s]

POD

DMD

5~10 [s]

POD

ここまで来て，ペアができていたのは，2つの渦ができていた所為だと分かった．双子渦しか立ち上がらないこの時間では，ペアがない．

DMD

10~15 [s]

POD

DMD

終わりに

5年以上前から名前は知っていたのだけれど．．．色々手を動かして初めて，DMDの中身を，少しだけ，きちんと，理解できた気がする．predictionとかの枠組みもできるという触れ込みだが，果たして？

周辺に，

spDMD: sparsity-promoting DMD
mrDMD: multi-resolution DMD
eDMD: extended DMD
jetDMD: jet DMD
piDMD: physics-informed DMD

など（Schmid2022）．