はじめに

職場でDMDの知識が求められる案件があったので、アウトプットとして投稿。
参考はData-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Controlの第七章。久々に英語の書籍を読むのが辛かった。
今回はexact DMDのアルゴリズムを理解して、公開されている二次元円柱周りの流れ場データに対して、実際にDMDを適用することをゴールとする。

概要

Dynamical Sytems=力学系というのは、時間経過に伴い、ある法則にしたがって変化する系、また或いは、初期状態が与えられた時、その後の全ての状態量の変化が決定される系のことを指す。

その概念は、かつてポアンカレが天体の運動を研究していた時に端を発する。彼は2つの天体による運動は比較的単純であることを発見したが、3つ目の天体が加わると運動が非常に複雑になり、その結果として予測が困難になることを明らかにした。この発見は、カオス理論の初期の形態と見なされている。その後、この理論は古典機械システム、電子回路、流体、気候学、金融、経済など多岐にわたる分野での数学的モデリングへと発展した。

取り分け近年では、支配方程式が未知である一方で、観測データが十分にあるため、解析的なアプローチからData-Drivenなアプローチへと移行している。

Dynamical Systems

本章では、近年のData-Driven型Dynamical Systemsへの導入の準備として、力学系の表記法を導入し、主要な動機と課題について述べる。

基本定義

初めに、Dynamical Systemsの数学的記述を導入する。
時間を$t,s\in T$、空間を$\boldsymbol{x}\in H$としたとき、

\begin{cases}
\text{恒等性：}&\boldsymbol{f}(\boldsymbol{x}, 0) = \boldsymbol{x}_0 \\\\
\text{結合法則:}&\boldsymbol{f}(\boldsymbol{f}(\boldsymbol{x}, t), s) = \boldsymbol{f}(\boldsymbol{x}, t+s)
\end{cases}

この条件を満足する接バンドル$\boldsymbol{f}$を含めた組み合わせを力学系と呼ぶ。参考書では、写像$\boldsymbol{f}$がリプシッツ連続であることを前提に話が進むので、上記2つの定義から以下のように一階微分方程式が容易に求められる。

\begin{align}
\boldsymbol{x}(t)&=\boldsymbol{f}(\boldsymbol{x}_0, t) \\
\Rightarrow \frac{\partial}{\partial t}\boldsymbol{x}(t)&=\frac{\partial}{\partial t}\boldsymbol{f}(\boldsymbol{x}_0, t) \\
&=\left.\frac{\partial\boldsymbol{f}}{\partial t}\right|_{(\boldsymbol{x}, t)}
\end{align}

参考書では、これをより簡潔に（というか物理の慣習的に?）$\frac{d}{dt}\boldsymbol{x}=\boldsymbol{f}(\boldsymbol{x}, t; \boldsymbol{\beta})$と記述されている。
実例として、ローレンツ系の方程式を考え、pytonで可視化してみる。

支配方程式

\begin{cases}
\frac{\partial x}{\partial t}=\sigma(y-x)\\
\frac{\partial y}{\partial t}=x(\rho -z)-y\\
\frac{\partial z}{\partial t}=xy-\beta z
\end{cases}

可視化コード（python）

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def lorenz(x_y_z, _, sigma, beta, rho):
    x, y, z = x_y_z
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

# Initial conditions
x0 = (0., 1., 20.)
eps = 0.01
t = np.arange(0, 50, eps)

# Lorenz attractor parameters
sigma = 10.0
beta = 2.
rho = 25.0

# Solving the Lorenz system
x_t = odeint(lorenz, x0, t, args=(sigma, beta, rho))

# Extracting individual components
x, y, z = x_t.T

# Creating a 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, linewidth=1)
plt.show()

可視化結果

Discrete-Time Systems

前章では連続力学系について紹介したが、本章は離散力学系について述べる。
離散力学系は、実験データやデジタルコントロールの分野ではより重要な系（マシンが扱えるから）となる。

基本定義

基本定義は以下の漸化式として表現できる。

\boldsymbol{x}_{k+1}=\boldsymbol{F}(\boldsymbol{x}_{k})

この離散力学系の基本定義は前述の連続力学系の基本定義から導出することが出来る。
十分に小さなステップ時間$\Delta t$を取ることで、多様体$\boldsymbol{x}_{k}=\boldsymbol{x}(k\Delta t)$となるので、

\begin{align}
&&{\text 連続力学系基本定義: }\frac{d}{dt}\boldsymbol{x}(t)&=\boldsymbol{f}(\boldsymbol{x}(t)) \\
&\Rightarrow&\boldsymbol{x}(t)&=\int_{0}^{t}\boldsymbol{f}(\boldsymbol{x}(\tau))d\tau \\
&\Rightarrow&\boldsymbol{x}(k\Delta t+\Delta t)&=\int_{0}^{k\Delta t+\Delta t}\boldsymbol{f}(\boldsymbol{x}(\tau))d\tau \\
&&&=\int_{0}^{k\Delta t}\boldsymbol{f}(\boldsymbol{x}(\tau))d\tau+\int_{k\Delta t}^{k\Delta t+\Delta t}\boldsymbol{f}(\boldsymbol{x}(\tau))d\tau \\
&&&= \boldsymbol{x}(k\Delta t)+\int_{k\Delta t}^{k\Delta t+\Delta t}\boldsymbol{f}(\boldsymbol{x}(\tau))d\tau \\
&i.e.&\boldsymbol{x}(t_o+t)=\boldsymbol{F}_t(\boldsymbol{x}(t_0))&=\boldsymbol{x}(t_0)+\int_{t_0}^{t_0+t}\boldsymbol{f}(\boldsymbol{x}(\tau))d\tau
\end{align}

この離散力学系は$\Delta t$を小さくすればする程、連続力学系へ近づき、計算機科学の分野では一般的な概念となっている。

線形力学系とスペクトル分解

線形力学系は閉形式の解を持つので、可能な限りそれとして扱うのが推奨される。

\begin{align}
\frac{d}{dt}\boldsymbol{x}(t)&=\boldsymbol{Ax}(t) \\
\Rightarrow \boldsymbol{x}(t_0+t)&=e^{\boldsymbol{A}t}\boldsymbol{x}(t_0)
\end{align}

ここで、行列指数関数$e^{\boldsymbol{A}}:=\sum^\infty_k \frac{1}{k!}A^k$である。
この力学系の特性は、行列$\boldsymbol{A}$の固有行列および固有値に依存する。行列$\boldsymbol{A}$のスペクトル分解（固有値分解）は、

\boldsymbol{AT}=\boldsymbol{T\Lambda}

と表現される。
ただし、行列$\boldsymbol{A}$は$n$個の相違なる固有値を持ち、行列$\boldsymbol{\Lambda}$は対角成分に$\lambda_j$を持つ対角行列、行列$\boldsymbol{T}$は固有値$\lambda_j$に付随する線形独立固有ベクトル$\xi_j$を列に持つ。（即ち、$\boldsymbol{T}=(\xi_1, \xi_2,...,\xi_n)$）
特に、上記の場合は行列$\boldsymbol{A}$が対角化可能であるので、$\boldsymbol{A}=\boldsymbol{T\Lambda T}^{-1}$と記述出来るので、

\begin{align}
\boldsymbol{x}(t_0+t)&=e^{\boldsymbol{A}t}\boldsymbol{x}(t_0) \\
&=e^{\boldsymbol{T\Lambda T}^{-1}t}\boldsymbol{x}(t_0) \\
&=\left\{\sum^\infty_{k=0} \frac{1}{k!}(\boldsymbol{T\Lambda T}^{-1}t)^k\right\}\boldsymbol{x}(t_0) \\
&=\left\{\sum^\infty_{k=0} \frac{1}{k!}\boldsymbol{T}(\boldsymbol{\Lambda}t)^k\boldsymbol{T}^{-1}\right\}\boldsymbol{x}(t_0) \\
&=\boldsymbol{T}e^{\boldsymbol{\Lambda}t}\boldsymbol{T}^{-1}\boldsymbol{x}(t_0)
\end{align}

重複する固有値を持つ際には、ジョルダン標準を用いることで表現可能だが、参考書第七章では解説されていないため、割愛する。
ところで、上記式に置いて両辺左から$\boldsymbol{T}^{-1}$を掛けると

\boldsymbol{T}^{-1}\boldsymbol{x}(t_0+t)=e^{\boldsymbol{\Lambda}t}\boldsymbol{T}^{-1}\boldsymbol{x}(t_0)

さらに、$\boldsymbol{z}=\boldsymbol{T}^{-1}\boldsymbol{x}$と置換して

\begin{align}
&&\boldsymbol{z}(t_0+t)&=e^{\boldsymbol{\Lambda}t}\boldsymbol{z}(t_0) \\
&\Rightarrow&\frac{d}{dt}\boldsymbol{z}&=\boldsymbol{\Lambda}\boldsymbol{z} \\
&\Rightarrow&\frac{d}{dt}z_j&=\lambda_j z_j\quad\because \boldsymbol{\Lambda}\text{ is a diagonal matrix}
\end{align}

このように、はじめに導入した線形力学系の式へと回帰できた。
即ち、ベクトル$\boldsymbol{x}$を固有空間へマッピングすることで、次元ごとの単純な力学系で表現することが可能である。

Dynamic Mode Decomposition (DMD)

Dynamic Mode Decomposition（以下、DMD）は、シュミットが流体力学分野で開発した、高次元データから時空間構造を抽出する次元削減手法である。そのアルゴリズムは特異値分解（SVD）や固有直交分解（POD）に基づくものだが、これらの手法が空間的な相関だけ基づく階層化を行うのに対して、DMDは時間情報を考慮に入れることが出来る。つまり、DMDアルゴリズムは、空間的次元削減のためのSVDと、時間的周波数同定のための高速フーリエ変換（FFT）のいいとこ取りアルゴリズムなのだ。しかし、そのトレードオフとして制約的になることに留意する必要がある。

特に線形システムの場合、DMDはクープマン作用素の近似となり、非線形力学系を線形演算により近似的に解析することが出来る。さらに、実験データや数値データに等しく有効であり、支配方程式に依存しないデータ駆動型の特性評価ができる。

この様な特性や、数値的な柔軟性もあり、DMDは近年では非常に人気が高まっていると言う。

DMDアルゴリズム概要

DMDアルゴリズムはいくつか提案されているが、ここではTuらによるexact DMDを紹介する。従来手法がデータの一様サンプリングを要求してくるのに対し、exact DMDは不均一サンプリングにも対応できる。

早速、数学的定式化を導入する。時間展開する系のスナップショットの組を以下のように表す。

\text{snapshot}:=\{\boldsymbol{x}(t_k), \boldsymbol{x}(t'_k)\}_{k=1}
^m

ここで、系の挙動を表現するのに十分に小さい$\Delta t$を用いて、漸化式$t'_k=t_k+\Delta t$が満足されるものとする。これらの列ベクトルを束ねて、次のような行列を定義する。

\begin{align}
\boldsymbol{X}=
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}(t_1) & \boldsymbol{x}(t_2) & \cdots & \boldsymbol{x}(t_m) \\
| & | & & |
\end{bmatrix}\\\\
\boldsymbol{X}'=
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}(t'_1) & \boldsymbol{x}(t'_2) & \cdots & \boldsymbol{x}(t'_m) \\
| & | & & |
\end{bmatrix}
\end{align}

さて、改めてDMDの目的は、以下に式を最も満足する線形演算$\boldsymbol{A}$の固有値、固有ベクトルを求めることである。固有値は複素数であり、その実部はモードの成長or減衰の率を表し、虚部はモードの振動周波数を示す。したがって、固有値を解析することで、システムが時間と共にどのように発展していくかを把握できる。

\boldsymbol{X}'\simeq\boldsymbol{AX}

シュミットによる元の定式化では一様サンプリングを仮定しており、その場合$\boldsymbol{x}_k=\boldsymbol{x}(k\Delta t)$と記述できるので、

\begin{align}
&&\boldsymbol{X}'&\simeq\boldsymbol{AX}\\
&\Rightarrow&
\left(
\boldsymbol{x}(t_1)\ \boldsymbol{x}(t_2)\ \cdots\ \boldsymbol{x}(t_m) \right)
&\simeq\boldsymbol{A}
\left(
\boldsymbol{x}(t_2)\ \boldsymbol{x}(t_3)\ \cdots\ \boldsymbol{x}(t_{m+1}) \right)\\
&\Rightarrow&\boldsymbol{x}_{k+1}&\simeq\boldsymbol{Ax}_k
\end{align}

即ち、最適な$\boldsymbol{A}$は以下のように与えられる。

\newcommand{\argmax}{\mathop{\rm arg~max}\limits}
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
\boldsymbol{A}=\argmin_{\boldsymbol{A}}\|\boldsymbol{X}'-\boldsymbol{AX}\|_F
=\boldsymbol{X}'\boldsymbol{X}^\dagger

ただし、$\|\cdot\|_F$はフロベニウスノルム、$^\dagger$は擬似逆行列を表す。

フロベニウスノルムについてはこちら

フロベニウスノルムは2つの行列がどのくらい違っているかを測るための尺度のひとつである。行列$\boldsymbol{A}\in\mathbb{C}^{n\times m}, n\geq m$を考えたとき、その特異値も$\sigma_k$と表現して

\|\boldsymbol{A}\|_F=\sqrt{\sum_i^m\sum_j^n\mid a_{ij}\mid^2}
=\sqrt{\mathrm{tr}(A^*A)}=\sqrt{\sum_{i=1}^m \sigma^2_i}

このノルムはフロベニウスノルムと呼ばれ、これを目的関数とした場合の最小化問題の解がSVD(what we call truncated SVD)である、という内容がEckart-Youngの定理である。

\newcommand{\argmax}{\mathop{\rm arg~max}\limits}
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
\argmin_{\tilde{\boldsymbol{X}}, \mathrm{s.t. rank}(\tilde{\boldsymbol{X}})=r}\|\boldsymbol{X}-\tilde{\boldsymbol{X}}\|_F
=\tilde{\boldsymbol{U}}\tilde{\boldsymbol{\Sigma}}\tilde{\boldsymbol{V}}^*

ここでは深掘りは行わないが、上記で吟味された行列$\boldsymbol{A}$が後述するクープマン作用素（実測定関数$\boldsymbol{g}(\boldsymbol{x})=\boldsymbol{x}$の場合）に似ていることに着目する。$\boldsymbol{A}$がクープマン作用素の線形測定有限部分空間の近似であることから、やはり行列$\boldsymbol{A}$についての固有値$\Lambda$と固有ベクトル行列$\boldsymbol{\Phi}$を求めることが目的となる。

\boldsymbol{A\Phi}=\boldsymbol{\Phi\Lambda}

exact DMDアルゴリズム

前節で言及した式のように、行列$\boldsymbol{A}$から固有値、固有ベクトルを求めることは理論的には単純であるが、リアルワールドデータのように高次元になってくると、スペクトル分解どころか、$\boldsymbol{A}$を正しく表現することすらままならないと言う。

DMDアルゴリズムは次元削減を利用して、$\boldsymbol{A}$を直接使った明示的な計算を必要とせずに、スペクトル分解を実現する。さらに、ある特定の条件下では、これらの縮小行列による近似が行列$\boldsymbol{A}$の固有ベクトルに完全一致する。(exact DMD)

以下より、その手順について解説する。

Step 1.

はじめに、系の状態$\boldsymbol{X}\in\mathbb{C}^{n\times m}$のSVDを行う。
部分ユニタリ行列$\tilde{\boldsymbol{U}}\in\mathbb{C}^{n\times r}$と$\tilde{\boldsymbol{V}}\in\mathbb{C}^{r\times m}$、非負実対角行列$\tilde{\boldsymbol{\Sigma}}\in\mathbb{C}^{r\times r}$を用いて、

\boldsymbol{X}\simeq\tilde{\boldsymbol{U}}\tilde{\boldsymbol{\Sigma}}\tilde{\boldsymbol{V}}^*

多くのケースで、$\boldsymbol{X}$はtall-skinnyな行列であり、階数$r$をどう選ぶかもDMDの過程でひとつ問題となる。（参考:Gavish and Donoho）

Step. 2

前節より、行列$\boldsymbol{A}$は擬似逆行列$\boldsymbol{X}$を用いて得ることができ、

\begin{align}
\boldsymbol{A}&=\boldsymbol{X}'\boldsymbol{X}^\dagger \\
&=\boldsymbol{X}'\left(\tilde{\boldsymbol{U}}\tilde{\boldsymbol{\Sigma}}\tilde{\boldsymbol{V}}^*\right)^{*}\\
&=\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\tilde{\boldsymbol{U}}^*
\end{align}

ここで、興味があるのは行列$\boldsymbol{A}$の固有値であるので、$\boldsymbol{X}$の左特異ベクトル縮約行列$\tilde{\boldsymbol{U}}$による行列$\boldsymbol{A}$の類似変換をすれば、

\begin{align}
\tilde{\boldsymbol{A}}&=\tilde{\boldsymbol{U}}^*\boldsymbol{A}\tilde{\boldsymbol{U}}\\
&=\tilde{\boldsymbol{U}}^*\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\tilde{\boldsymbol{U}}^*\right)\tilde{\boldsymbol{U}}\\
&=\tilde{\boldsymbol{U}}^*\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\quad\because\boldsymbol{U}^*\boldsymbol{U}=\boldsymbol{E}
\end{align}

このように、同じ$\boldsymbol{A}$と同じ固有値を持つ縮約行列$\tilde{\boldsymbol{A}}$をより単純な行列の積から得ることができる。

類似変換の固有値保存についてはこちら

部分ユニタリ行列$\boldsymbol{U}\in\mathbb{C}^{n\times r}$による類似変換の永年方程式の左辺を展開する。

\begin{align}
|\boldsymbol{U}^*\boldsymbol{A}\boldsymbol{U}-\lambda\boldsymbol{E}|
&=|\boldsymbol{U}^*\boldsymbol{A}\boldsymbol{U}-\lambda\boldsymbol{U}^*\boldsymbol{U}|\\
&=|\boldsymbol{U}^*\left(\boldsymbol{A}-\lambda\boldsymbol{E}\right)\boldsymbol{U}|\\
&=|\boldsymbol{A}-\lambda\boldsymbol{E}|=0
\end{align}

よって、固有値が保存されることがわかる。ただし、固有ベクトルの保存については保証されない。

行列$\boldsymbol{U}$による縮約状態へのマッピング$\tilde{\boldsymbol{x}}=\tilde{\boldsymbol{U}}^*\boldsymbol{x}$を実施することで、縮約行列$\boldsymbol{A}$による線形ダイナミクス$\tilde{\boldsymbol{x}}_{k+1}=\tilde{\boldsymbol{A}}\tilde{\boldsymbol{x}}_k$に置き換えることができる。（次元削減）

Step. 3

前Step.で求めた縮約行列$\tilde{\boldsymbol{A}}$のスペクトル分解を計算する。

\tilde{\boldsymbol{A}}\boldsymbol{W}=\boldsymbol{W\Lambda}

対角行列$\boldsymbol{\Lambda}$の対角成分はDMDの固有値であり、元の行列$\boldsymbol{A}$のそれと一致する。固有ベクトル行列$\boldsymbol{W}$は$\tilde{\boldsymbol{A}}$の固有ベクトルからなる行列で、対角化する座標変換を行う。

Step. 4

さて、ここで漸くDMDの固有ベクトル行列$\boldsymbol{\Phi}$が導出できるようになる

\begin{align}
&&\tilde{\boldsymbol{A}}\boldsymbol{W}
&=\boldsymbol{W\Lambda}\\
&\Rightarrow&\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\right)\tilde{\boldsymbol{A}}\boldsymbol{W}
&=\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\right)\boldsymbol{W\Lambda}\\
&\Rightarrow&\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\right)\left(\tilde{\boldsymbol{U}}^*\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\right)\boldsymbol{W}
&=\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\right)\boldsymbol{W}\boldsymbol{\Lambda}\\
&\Rightarrow&\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\tilde{\boldsymbol{U}}^*\right)\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{W}\right)
&=\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{W}\right)\boldsymbol{\Lambda}\\
&\Rightarrow&\boldsymbol{A}\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{W}\right)
&=\left(\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{W}\right)\boldsymbol{\Lambda}\\
&\Rightarrow&\boldsymbol{\Phi}&=\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{W}
\end{align}

シュミットの元の論文中では、$\boldsymbol{\Phi}=\tilde{\boldsymbol{U}}\boldsymbol{W}$とされており、projected modesと呼ばれている。先ほど導出した固有ベクトル行列$\boldsymbol{\Phi}$が$\boldsymbol{X}'$の列空間に存在するのに対し、projected modesは$\boldsymbol{X}$の列空間に存在することになる。線形変換$\boldsymbol{A}$の定義は$\boldsymbol{A}:=\boldsymbol{X}'\boldsymbol{X}^\dagger$であることから、本来であれば$\boldsymbol{\Phi}$は$\boldsymbol{X}'$の列空間にあるべきだが、実際の$\boldsymbol{X}$と$\boldsymbol{X}'$の列空間はほぼ同じ（直感的で良い）なので、あまり大きな問題にはならない。

従来のDMDと本書籍が解説に注力しているexact DMDの違いはまさにここで、数学的な正確性と計算コストのトレードオフで前者を取ったのが、exact DMDである。

固有値$\lambda=0$に対するDMDモードが$\boldsymbol{\phi}=\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{w}=0$の場合、$\boldsymbol{\phi}=\tilde{\boldsymbol{U}}\boldsymbol{w}$を用いると良い。

歴史的背景

これは余談として記されているのだが、元々$\boldsymbol{X}$と$\boldsymbol{X}'$は以下のように、連続したスナップショットをズラした形で記述されていた。

\begin{align}
\boldsymbol{X}&=
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}_1 & \boldsymbol{x}_2 & \cdots & \boldsymbol{x}_m \\
| & | & & |
\end{bmatrix}\\\\
\boldsymbol{X}'&=
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}_2 & \boldsymbol{x}_3 & \cdots & \boldsymbol{x}_{m+1} \\
| & | & & |
\end{bmatrix}
\end{align}

つまり、$\boldsymbol{X}$は行列$\boldsymbol{A}$の反復で、

\boldsymbol{X}\simeq
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}_1 & \boldsymbol{A}\boldsymbol{x}_1 & \cdots & \boldsymbol{A}^{m-1}\boldsymbol{x}_1 \\
| & | & & |
\end{bmatrix}

と近似で表現できる。
ここで、各列がクリロフ部分空間（c.f.Wikipedia）に属することに注目する。行列$\boldsymbol{X}'$は$\boldsymbol{X}$をシフト演算した形で表される。

\begin{align}
\boldsymbol{X}'&=\boldsymbol{XS} \\
\boldsymbol{S}&:=
\begin{bmatrix}
0 & 0 & \cdots & 0 & a_1 \\
1 & 0 & \cdots & 0 & a_2 \\
0 & 1 & \cdots & 0 & a_3 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots & 1 & a_m \\
\end{bmatrix}
\end{align}

行列$\boldsymbol{S}$の形を見てみると、左$m-1$列は列のシフト、一番右の列は残差を最小化する係数が与えられている。
このように、当初のDMDアルゴリズムはアーノルディ法（固有値計算反復アルゴリズム）とよく似ていることが分かる。実際に、$\boldsymbol{X}=\left(\boldsymbol{x}_1\ \boldsymbol{x}_2\ \cdots\ \boldsymbol{x}_m\right)$はクリロフ部分空間$\mathcal{K}_m(\boldsymbol{A},\boldsymbol{x}_1)$を構成する列ベクトルからなる行列であり、これを列正規化したものを$\boldsymbol{Q}$とすると、

\begin{align}
&\boldsymbol{Q}=\left(\frac{\boldsymbol{x}_1}{\|\boldsymbol{x}_1\|}\ \frac{\boldsymbol{x}_2}{\|\boldsymbol{x}_2\|}\  \cdots\ \frac{\boldsymbol{x}_m}{\|\boldsymbol{x}_m\|}\right) \\
\Rightarrow &\boldsymbol{S}\simeq\boldsymbol{Q}^*\boldsymbol{A}\boldsymbol{Q} \\
\Rightarrow &\boldsymbol{Q}\boldsymbol{S}\simeq\boldsymbol{A}\boldsymbol{Q}
\end{align}

ここで、$\boldsymbol{S}$の固有値を$\lambda$、固有ベクトルを$\boldsymbol{v}$、誤差項を$\boldsymbol{D}$として、

\begin{align}
&\boldsymbol{Q}\boldsymbol{S}\simeq\boldsymbol{A}\boldsymbol{Q} \\
\Rightarrow&\boldsymbol{Q}\boldsymbol{S}+\boldsymbol{D}=\boldsymbol{A}\boldsymbol{Q} \\
\Rightarrow&\boldsymbol{Q}\boldsymbol{S}\boldsymbol{v}+\boldsymbol{D}\boldsymbol{v}=\boldsymbol{A}\boldsymbol{Q}\boldsymbol{v} \\
\Rightarrow&\lambda\boldsymbol{Q}\boldsymbol{v}+\boldsymbol{D}\boldsymbol{v}=\boldsymbol{A}\boldsymbol{Q}\boldsymbol{v} \\
\Rightarrow&\boldsymbol{D}\boldsymbol{v}=\left(\boldsymbol{A}-\lambda\boldsymbol{E}\right)\boldsymbol{Q}\boldsymbol{v} \\
\end{align}

誤差$\boldsymbol{D}\rightarrow\boldsymbol{O}$の時、$\boldsymbol{A}-\lambda\boldsymbol{E}\rightarrow\boldsymbol{O}$なので、$\boldsymbol{S}$と$\boldsymbol{A}$の固有値は共有されることになる。したがって、行列$\boldsymbol{S}$のスペクトル分解をすることでDMDが実現できるのだが、先ほど紹介したexact DMDアルゴリズムの方が数値計算的に優れている。

特異値分解とスペクトル展開

DMDの重要なポイントは、データドリブンなスペクトル分解という観点からシステム状態を展開できることだ。
$\boldsymbol{\phi}_j$をDMDモード（行列$\boldsymbol{A}$の固有ベクトル）、$\lambda_j$をDMD固有値（行列$\boldsymbol{A}$の固有値）として、

\begin{align}
\boldsymbol{x}_k&=\boldsymbol{A}^{k-1}\boldsymbol{x}_1 \\
&=\boldsymbol{\Phi\Lambda}^{k-1}\boldsymbol{\Phi}^\dagger\boldsymbol{x}_1
\end{align}

ここで、モード振幅$b_j$を$\left(b_0\ b_1\cdots b_r\right)=\boldsymbol{\Phi}^\dagger\boldsymbol{x}_1$を満足するものと置けば、

\begin{align}
\boldsymbol{x}_k&=\boldsymbol{\Phi}\boldsymbol{\Lambda}^{k-1}\boldsymbol{b}　 \\
&=
\begin{bmatrix}
| & & | \\
\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_r \\
| & & |
\end{bmatrix}
\begin{bmatrix}
\lambda_1^{k-1} & & \\
& \ddots & \\
& & \lambda_r^{k-1}
\end{bmatrix}
\begin{bmatrix}
b_1 \\
\vdots \\
b_r
\end{bmatrix} \\
&=\sum_{j=1}^r \boldsymbol{\phi}_j\lambda_j^{k-1}b_j
\end{align}

参考書では、

\boldsymbol{x}_k=\boldsymbol{\Phi}\boldsymbol{\Lambda}^{k-1}\boldsymbol{b}=
\begin{bmatrix}
| & & | \\
\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_r \\
| & & |
\end{bmatrix}
\begin{bmatrix}
\lambda_1 & & \\
& \ddots & \\
& & \lambda_r
\end{bmatrix}
\begin{bmatrix}
b_1 \\
\vdots \\
b_r
\end{bmatrix}

とされているが、固有行列の指数が抜けている気がする。

と各モード$\boldsymbol{x}_k$を特異値分解することができる。
さらに、このDMDスペクトル展開は以下のように式変形可能

\begin{align}
\boldsymbol{x}_k&=
\begin{bmatrix}
| & & | \\
\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_r \\
| & & |
\end{bmatrix}
\begin{bmatrix}
b_1 & & \\
& \ddots & \\
& & b_r
\end{bmatrix}
\begin{bmatrix}
\lambda_1^{k-1} \\
\vdots \\
\lambda_r^{k-1}
\end{bmatrix} \\
\Rightarrow
\boldsymbol{X}&=
\begin{bmatrix}
| & & | \\
\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_r \\
| & & |
\end{bmatrix}
\begin{bmatrix}
b_1 & & \\
& \ddots & \\
& & b_r
\end{bmatrix}
\begin{bmatrix}
\lambda_1 & \cdots & \lambda_1^{m-1} \\
\vdots & \ddots & \vdots \\
\lambda_r & \cdots & \lambda_r^{m-1}
\end{bmatrix}
\end{align}

モード振幅$\boldsymbol{b}$は最初のスナップショット$\boldsymbol{x}_1$を使ってその混合を決定するのだが、$\boldsymbol{b}=\boldsymbol{\Phi}^\dagger\boldsymbol{x}_1$という簡単な定義式を用いても計算コストがかなり高くなってしまう。
そこで、projected modes同様に縮約モデルを用いた計算コスト軽減を図る。

\begin{align}
&&\boldsymbol{x}_1&=\boldsymbol{\Phi b} \\
&\Rightarrow&\tilde{\boldsymbol{U}}\tilde{\boldsymbol{x}}_1&=\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{Wb} \\
&\Rightarrow&\tilde{\boldsymbol{x}}_1&=\tilde{\boldsymbol{U}}^*\boldsymbol{X}'\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Sigma}}^{-1}\boldsymbol{Wb} \\
&\Rightarrow&\tilde{\boldsymbol{x}}_1&=\tilde{\boldsymbol{A}}\boldsymbol{Wb} \\
&\Rightarrow&\tilde{\boldsymbol{x}}_1&=\boldsymbol{W\Lambda b} \\
&\Rightarrow&\boldsymbol{b}&=\left(\boldsymbol{W\Lambda}^{-1}\right)\tilde{\boldsymbol{x}}_1 \\
\end{align}

元の固有ベクトル行列$\boldsymbol{\Phi}\in\mathbb{R}^{n\times r}$と比較して、行列$\boldsymbol{W}$と$\boldsymbol{\Lambda}$はいずれも$\mathbb{R}^{r\times r}$に属しており、軽量化されていることが分かる。

上記の特異値展開は、連続固有値$\omega=\log{(\lambda)/\Delta t}$を導入することによって、連続時間でも記述できる。

\boldsymbol{x}(t)=\sum_{j=1}^r\boldsymbol{\phi}_je^{\omega_j t}b_j=\boldsymbol{\Phi}e^{\boldsymbol{\Omega}t}\boldsymbol{b}

ここで、$\boldsymbol{\Omega}$は連続系の固有値$\omega_j$を対角成分に持つ対角行列であり、データ行列$\boldsymbol{X}$はVandermonde行列$\boldsymbol{T}(\boldsymbol{\omega})$を用いて以下のように表現できる。

\boldsymbol{X}\simeq
\begin{bmatrix}
| & & | \\
\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_r \\
| & & |
\end{bmatrix}
\begin{bmatrix}
b_1 & & \\
& \ddots & \\
& & b_r
\end{bmatrix}
\begin{bmatrix}
e^{\omega_1t_1} & \cdots & e^{\omega_1t_m} \\
\vdots & \ddots & \vdots \\
e^{\omega_rt_1} & \cdots & e^{\omega_rt_m}
\end{bmatrix}
=\boldsymbol{\Phi}\text{diag}(\boldsymbol{b})\boldsymbol{T}(\boldsymbol{\omega})

コーディング例

本章では、DMDのプログラム実装を行う。
シミュレーションデータはIBPMで作成された、二次元円柱周りの流れ場データを用いる。データはこちらから取得可能。

環境

言語: python3.12
OS: macOS Sonoma14.5
RAM: 64GB

データ読み込み

Jupyterで書いたので、セル別に区切って掲載する。

import time

import matplotlib.animation as animation
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
from scipy import io
from tqdm import tqdm

vortall_mat = io.loadmat('../DMD/FLUIDS/CYLINDER_ALL.mat')
X = vortall_mat['VORTALL']
time_steps = X.shape[1]
width = 199
height = 449

FLUIDSというディレクトリ傘下にCYLINDER_ALL.matというMatlab配列のデータがあるので、それを読み取りに行く。以下は必要ないのだが、vortall_mat.keys()でキーを確認すると、'UALL', 'UEXTRA', 'VALL', 'VEXTRA', 'VORTALL', 'VORTEXTRA'があることが分かる。末尾がEXTRAのものはデータ数が1でALLは151あるので、（推測だが）EXTRAは時系列的に一単位時間未来のデータなのだろう。今回は'VORTALL'を使う。

シミュレーションデータの可視化

DMDをする前に、まずは元のシミュレーションデータの可視化をしたくなる。Matplotlibには様々なカラーマップが用意されているが、離散で配色がグラデーションとなっているものがないので、カスタムで作成する。

cmap = mcolors.ListedColormap([
    '#67001f', '#7f1825', '#98312b', '#b2182b', '#ca3a4b',
    '#d6604d', '#e78a70', '#f4a582', '#fddbc7', '#fbf0ef',
    '#f7f7f7', '#d1e5f0', '#b1cfe9', '#92c5de', '#6d9bc4',
    '#4393c3', '#2c7bae', '#2166ac', '#174a88', '#053061'
])

このカラーマップを用いてアニメーションを作成

### Animation of the original data
fig, ax = plt.subplots(figsize=(2.5, 4))
im = ax.imshow(
    X[:,0].reshape(height, width),
    cmap=cmap,
    norm=mcolors.Normalize(vmin=-3., vmax=3.),
    interpolation='nearest'
)
plt.colorbar(im, ax=ax)

### Update function for animation
def update(frame):
    im.set_array(X[:,frame].reshape(height, width))
    return [im]

ani = animation.FuncAnimation(
    fig,
    update,
    frames=range(time_steps),
    blit=True,
    interval=50
)
ani.save('flow_field_animation.gif', writer='imagemagick')

すると、以下のような.GIFが作成される。

DMDアルゴリズムの実装

先述したexact DMDの記述に倣って関数を作成する。
コメントアウトがある段落はアルゴリズムの各stepに対応している。

def DMD(X, Xprime, rank, exact = True):
    ### Compute the SVD of X
    U, S, Vh = np.linalg.svd(X, full_matrices=False)
    Ur = U[:,:rank]
    Sr = np.diag(S[:rank])
    Vhr = Vh[:rank, :]

    ### Obtain Atilde by computing the pseudo-inverse of X
    Atilde = np.linalg.solve(Sr, (Ur.T @ Xprime @ Vhr.T))

    ### The Spectral decomposition of Atilde
    Lambda, W = np.linalg.eig(Atilde)
    Lambda = np.diag(Lambda)

    ### Reconstructed the high-demensional DMD modes
    Phi = Xprime @ np.linalg.solve(Sr.T, Vhr).T @ W if exact else Ur @ W
    x1tilde = Sr @ Vhr[:,0]
    b = np.linalg.solve(W @ Lambda, x1tilde)
    return Phi, Lambda, b

そして、この関数を呼び出した結果を用いて（今回は$rank=21$）①DMDモード、②振動周波数を並べたアニメーションを作成

### Animation of the DMD modes (Diagnostics)
vmax = 10**-2
Phi, Lambda, b = DMD(X[:,:-1], X[:,1:], 21)

### Plot DMD modes
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
im1 = ax1.imshow(
    Phi[:,0].real.reshape(height, width),
    cmap=cmap,
    norm=mcolors.Normalize(vmin=-vmax, vmax=vmax),
    interpolation='nearest'
)
plt.colorbar(im, ax=ax1)

### ax2 customize
ax2.set_xlim(0, time_steps)
ax2.set_ylim(-120, 120)

### Plot oscillation frequencies
omega = np.log(np.diag(Lambda))
time_dynamics = np.zeros((len(b), time_steps), dtype='complex')
for t in range(time_steps):
    time_dynamics[:,t] = np.diag(b) @ np.exp(omega * t).T
line_real, = ax2.plot(np.arange(time_steps), time_dynamics[0].real, label='Real')
line_imag, = ax2.plot(np.arange(time_steps), time_dynamics[0].imag, label='Imaginary')

### Update function
def update(frame):
    im1.set_array(Phi[:,frame].real.reshape(height, width))
    im1.set_norm(mcolors.Normalize(vmin=-vmax, vmax=vmax))
    line_real.set_ydata(time_dynamics[frame].real)
    line_imag.set_ydata(time_dynamics[frame].imag)
    ax2.legend(loc='upper right')
    fig.suptitle(f"mode={frame+1}")
    return [im1, line_real, line_imag]

ani = animation.FuncAnimation(
    fig,
    update,
    frames=range(21),
    blit=True,
    interval=400
)
ani.save('mode_real_animation.gif', writer='imagemagick')

実行結果は以下のようになる。

実部を青線、虚部を橙線で示しているが、改めて述べると、実部はモードの成長率または減衰率、虚部はモードの振動周波数を表す。 故に、この結果を見ればモード18,19がシステム内でとりわけ大きなエネルギーを持っていると解釈できる。

システムの再構築

今度は逆に、DMDの結果を用いて系自体を再構築してみる。先述してはいるが、数学的に以下の操作を行うことになる。

\boldsymbol{X}\simeq
\begin{bmatrix}
| & & | \\
\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_r \\
| & & |
\end{bmatrix}
\begin{bmatrix}
b_1 & & \\
& \ddots & \\
& & b_r
\end{bmatrix}
\begin{bmatrix}
e^{\omega_1t_1} & \cdots & e^{\omega_1t_m} \\
\vdots & \ddots & \vdots \\
e^{\omega_rt_1} & \cdots & e^{\omega_rt_m}
\end{bmatrix}
=\boldsymbol{\Phi}\text{diag}(\boldsymbol{b})\boldsymbol{T}(\boldsymbol{\omega})

ひとまず、時間$t=0$についてシステムの復元をしてみる。

### Reconstructed X(t=0)
Xr = Phi @ np.diag(b) @ np.exp(omega).T
fig, ax = plt.subplots(figsize=(2.5, 4))
im = ax.imshow(
    Xr.real.reshape(height, width),
    cmap=cmap,
    norm=mcolors.Normalize(vmin=-3., vmax=3.),
    interpolation='nearest'
)
plt.colorbar(im, ax=ax)

すると以下のような結果を得る。

しかしながら、これだけ見てもどのくらい再現できているか定性的・定量的に判断できないため、元のシミュレーションとの絶対誤差を取ってみて、差分を視覚化してみる。

### Calculate the Frobenius norm of the difference from the original data.
im = ax.imshow(
    np.abs(Xr.real-X[:,0]).reshape(height, width),
    cmap=cmap,
    norm=mcolors.Normalize(vmin=-3., vmax=3.),
    interpolation='nearest'
)
plt.colorbar(im, ax=ax)

カラーマップを見比べてみると、視覚的には20%ほど数値誤差が出ていそうだということがわかる。次に定量的な指標としてここの記事でも述べているが、階数$rank$の近似率は以下のように表せられる。

\text{appx. ratio}:=\frac{\|\boldsymbol{X}-\boldsymbol{\Phi}\text{diag}(\boldsymbol{b})\boldsymbol{T}(\boldsymbol{\omega})\|_F}{\|\boldsymbol{X}\|_F}

num = np.linalg.norm((Xr.real-X[:,0]).reshape(height, width), 'fro')
den = np.linalg.norm(X[:,0].reshape(height, width), 'fro')
num*100/den

この結果として14.902を得る。故に、定量的にも定性的にもまずまずといった結果を得られた。

exact DMD VS standard DMD

先の解説では、exact DMDは計算コストを犠牲に数学的正しさを優先した手法だと述べた。しかしながら実は、上で実装したDMD関数でexact=Falseとして通常のDMDを計算すると、誤差が0.075を得る。
おそらく、実際の計算では行列をいくつも掛ける方が累積誤差が多くなってしまうのだろう。もちろん、計算コストはexact DMDの方が掛かるはずであるが、下記コードで確認してみる。

### Measurement of execution time
def measure_execution_time(func, args=(), kwargs={}, num_trials=10):
    execution_times = []
    
    for _ in tqdm(range(num_trials)):
        start_time = time.time()
        func(*args, **kwargs)
        end_time = time.time()
        
        execution_time = end_time - start_time
        execution_times.append(execution_time)
    
    average_time = np.mean(execution_times)
    variance_time = np.var(execution_times)
    
    return average_time, variance_time

exact_ave, exact_var = measure_execution_time(
    DMD,
    args=(X[:,:-1],X[:,1:], 21),
    kwargs={'exact':True},
    num_trials=30
)
standard_ave, standard_var = measure_execution_time(
    DMD,
    args=(X[:,:-1],X[:,1:], 21),
    kwargs={'exact':False},
    num_trials=30
)
print(f'exact DMD 平均実行時間:{exact_ave:.3}sec 分散:{exact_var:.3}')
print(f'standard DMD 平均実行時間:{standard_ave:.3}sec 分散:{standard_var:.3}')

これの実行結果は以下の通り。

100%|██████████| 30/30 [00:29<00:00,  1.01it/s]
100%|██████████| 30/30 [00:26<00:00,  1.13it/s]
exact DMD 平均実行時間:0.988sec 分散:0.0137
standard DMD 平均実行時間:0.881sec 分散:0.00211

これらの結果から、実務では精度、計算コストの観点からstandard DMDを用いるのが良さそうである。
机上の数学と現実の計算機科学の違いを目の当たりにした。

固有値のプロット

最後の検証として、固有値$\lambda_j$の複素平面へのプロットを行う。

### Plot eigenvalues on complex plane
re = Lambda.real
im = Lambda.imag
plt.figure(figsize=(4,4))
plt.scatter(re, im, color='blue', marker='o')

### Add an unit circle
theta = np.linspace(0, 2*np.pi, 300)
x = np.cos(theta)
y = np.sin(theta)
plt.plot(x, y, 'k-', linewidth=0.5)

plt.gca().set_aspect('equal', adjustable='box')
plt.xticks(np.arange(-1, 1.5, 0.5))
plt.yticks(np.arange(-1, 1.5, 0.5))
plt.xlabel('Re(λ)')
plt.ylabel('Im(λ)')
plt.show()

結果は以下の通り。

綺麗に単位円上に固有値が並んでいることが分かる。
DMDは実験データ、数値データどちらにも適用可能であるが、ノイズの多いデータに対しては敏感であることが知られている。そこで、元のデータに$10^{-1}$程度のガウシアンノイズを加えて、再度プロットしてみる。

### Added Gaussian noise
noise = np.random.normal(0, 10**-1, X.shape)
X_gn = X + noise
Phi_gn, Lambda_gn, b_gn = DMD(X_gn[:,:-1], X_gn[:,1:], 21)

このように、ノイズが加えられたデータのDMD固有値は不規則に配置される。

おわりに

今回はexact DMDのアルゴリズムを理解して、二次元円柱周りの流れ場データを使って実際にDMDを適用させてみた。参考書にはDMDの派生や近似手法などの紹介もあったので、実際にそれらの手法も試してみたいと思う。