更新履歴

2025/10/18: 初版

はじめに

以前投稿した動的モード分解の中で、数学的に正確な離散時間Koopman作用素の近似となるExact DMDを紹介した。今回は、Exact DMD以外の主要な派生手法と用途をいくつか紹介する。

目次

Optimized DMD
Sparsity-promoting DMD
Randomized DMD

前提

各手法を解説する前に、前提となる部分を先に記述する。

表記と命名法

文献によってはエルミート転置を $^*$ などと表記したりするが、本稿では単純な共役とエルミート転置を区別するために、表記を以下の様に統一する：

$\boldsymbol{A}^\top$：転置
```
A.T
```
$\overline{\boldsymbol{A}}$：複素共役行列
```
A.conj()
```
$\boldsymbol{A}^*$：エルミート転置
```
A.conj().T
```
$\boldsymbol{A}^\dagger$：擬似逆行列
```
np.linalg.pinv(A)
```

DMDの数学的定式化

時間展開する系 $\boldsymbol{x}_k$ のスナップショットの組を以下のように表す：

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

多くの場合（例えばシミュレーションデータなど）では、時刻は等間隔に並んでおり、 $t'_k=t_k+\Delta t$ が成立する。故に、スナップショット行列は次のように構成できる：

\begin{align}
\boldsymbol{X}&=
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}_1 & \boldsymbol{x}_2 & \cdots & \boldsymbol{x}_{\tau-1} \\
| & | & & |
\end{bmatrix}
\tag{2}
\\
\boldsymbol{X}'&=
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}_2 & \boldsymbol{x}_3 & \cdots & \boldsymbol{x}_\tau \\
| & | & & |
\end{bmatrix}
\tag{3}
\end{align}

このとき、以下のような線形演算子 $\boldsymbol{A}$ を考える：

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

これは観測関数 $g(\boldsymbol{x})=\boldsymbol{x}$ の時の離散時間クープマン作用素の近似とみなせる。
$\boldsymbol{A}$ をスペクトル分解することで、力学系の特徴的な構造（モードや発展速度）を抽出できる。その固有値を $\lambda_j$、対応する固有ベクトルを $\boldsymbol{\phi}_j$ とすると：

\boldsymbol{A}\boldsymbol{\phi}_i=\lambda_i\boldsymbol{\phi}_i\tag{5}

このようにして得られる固有モードと固有値は、システムの振動数や成長率といった動的特徴を反映しており、時間発展は以下のように表現される：

\boldsymbol{x}_k=\sum_{j=1}^r \boldsymbol{\phi}_j\lambda_j^{k-1}b_j\tag{6}

ここで $b_j$ は初期状態に応じたモード振幅である。

また、本稿では行列演算が頻出するが、The Matrix Cookbookなどを参考にした。

データ準備

前回同様にシミュレーションデータにはIBPMで作成された、二次元円柱周りの流れ場データを用いる。データはこちらから取得可能。

環境

言語: python3.13
OS: macOS Sequoia 15.4
SoC: Apple M2 Max
RAM: 64GB

データ読み込み

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

import matplotlib.animation as animation
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
from scipy import io
from scipy.optimize import minimize
from scipy.signal import periodogram, find_peaks
from tqdm import tqdm

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

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

可視化のためにカスタムのカラーマップと描画関数を用意する。

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

# Create and return an animation of the time-varying field.
def animate_field(
    X, height, width, time_steps, cmap=cmap, vmin=-3.0, vmax=3.0, interval=50
):
    fig, ax = plt.subplots(figsize=(2.5, 4))
    im = ax.imshow(
        X[:, 0].reshape(height, width),
        cmap=cmap,
        norm=mcolors.Normalize(vmin=vmin, vmax=vmax),
        interpolation="nearest",
    )
    plt.colorbar(im, ax=ax)

    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=interval
    )
    return ani

これを用いて次の様にシミュレーションデータを可視化&保存できるようになる。

ani = animate_field(X, height, width, time_steps)
ani.save('cylinder_vorticity.gif', writer='imagemagick', fps=20)

[出力結果]

Optimized DMD (K.K.Chen et al., 2011)

Optimized DMDは非線形最小二乗最適化を用いて、Standard DMDより高精度なモード行列、固有値、振幅をデータ駆動で推定する手法である。

手法と数式モデル

Optimized DMDの基本的な考え方は次の通り。まず、固有値 $\lambda_j$ とモードベクトル $\boldsymbol{\psi}_j$、残差 $\boldsymbol{\varepsilon}_j$ を用いて、状態 $\lbrace\boldsymbol{x}_k\rbrace_{k=1}^\tau$ を以下のように表現する：

\boldsymbol{x}_k\simeq\sum_{j=1}^r\lambda_j^{k-1}\boldsymbol{\psi}_j+\boldsymbol{\varepsilon}_k\tag{7}

そして、次式を最小化する：

\Gamma:=\sum_{k=1}^{\tau} \|\boldsymbol{\varepsilon}_k\|_2^2\tag{8}

つまり、各時刻 $\boldsymbol{x}_k$ をモードの線型結合+残差と表現し、残差が小さくなるように固有値 $\lambda_j$ を決めるという考え方である。
各変数を行列に書き換えると次のようになる：

\begin{align}
\boldsymbol{\Psi}:&=
\begin{bmatrix}
    \boldsymbol{\psi}_1 & \cdots & \boldsymbol{\psi}_r
\end{bmatrix}\tag{9}
\\
\boldsymbol{V}_{\rm and}:&=
\begin{bmatrix}
    1 & \lambda_1 & \lambda_1^2 & \cdots & \lambda_1^{\tau-1} \\
    1 & \lambda_2 & \lambda_2^2 & \cdots & \lambda_2^{\tau-1} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    1 & \lambda_r & \lambda_r^2 & \cdots & \lambda_r^{\tau-1} \\
\end{bmatrix}\tag{10}
\\
\boldsymbol{E}:&=
\begin{bmatrix}
    \boldsymbol{\varepsilon}_1 & \cdots & \boldsymbol{\varepsilon}_r
\end{bmatrix}\tag{11}
\end{align}

これを用いれば、式 $(7)$ と式 $(8)$ は下記のように書き表せる：

\begin{align}
&
\begin{cases}
\boldsymbol{X}\simeq\boldsymbol{\Psi}\boldsymbol{V}_{\rm and}+\boldsymbol{E} \\
\min_{\lambda_j}\|\boldsymbol{E}\|_F^2
\end{cases}
\tag{12}
\\
\Rightarrow
&
\min_{\lambda_j}
\|
    \boldsymbol{X}-\boldsymbol{\Psi}\boldsymbol{V}_{\rm and}
\|_F^2
\tag{13}
\end{align}

Eckart-Youngの定理から、これを最小化する $\boldsymbol{\Psi}$ は：

\boldsymbol{\Psi}=\boldsymbol{XV}^\dagger_{\rm and}\tag{14}

そうすると、残差行列 $\boldsymbol{E}$ は：

\boldsymbol{E}=\boldsymbol{X}-\boldsymbol{XV}^\dagger_{\rm and}\boldsymbol{V}_{\rm and}\tag{15}

元々の最小化したいコスト関数 $\Gamma$ は次のように表せられる：

\begin{align}
\Gamma
=
\|\boldsymbol{X}-\boldsymbol{XV}^\dagger_{\rm and}\boldsymbol{V}_{\rm and}\|_F^2
=\|\boldsymbol{X}\|_F^2-\|\boldsymbol{XV}^\dagger_{\rm and}\boldsymbol{V}_{\rm and}\|_F^2\tag{16}
\end{align}

式 $(16)$ の変形について、フロベニウスノルムの分配法則は一般に成立しないが、1. $\boldsymbol{V}_{\rm and}$ はshort-fatな行列、つまり $m>r$ であること。2.基本的に固有値 $\lambda_j$ は重複しないこと、を考えれば、$\boldsymbol{V}_{\rm and}$ はフル行ランクとなり：

\boldsymbol{V}_{\rm and}^\dagger
=\boldsymbol{V}_{\rm and}^*
\left(
    \boldsymbol{V}_{\rm and}
    \boldsymbol{V}_{\rm and}^*
\right)^{-1}
\tag{17}

が成り立つ。このとき、行列 $\boldsymbol{P}=\boldsymbol{V}^\dagger_{\rm and}\boldsymbol{V}_{\rm and}\in\mathbb{C}^{\tau\times\tau}$ は $\boldsymbol{V}_{\rm and}$ の行空間への正射影行列となり、$\boldsymbol{XP}$ と $\boldsymbol{X(I-P)}$ が直行するので、式 $(17)$ のフロベニウスノルムの分配法則が成立する。

とは書いたものの、実際には $\boldsymbol{S}:=\boldsymbol{X}^*\boldsymbol{X}$ の計算が再利用できることから、トレースを用いた損失が計算では用いられる：

\begin{align}
\Gamma
&=
\|\boldsymbol{X}-\boldsymbol{XV}^\dagger_{\rm and}\boldsymbol{V}_{\rm and}\|_F^2
\tag{18}
\\
&=\|\boldsymbol{X}
(
    \boldsymbol{I}
    -\boldsymbol{V}^\dagger_{\rm and}\boldsymbol{V}_{\rm and}
)\|_F^2
\tag{19}
\\
&=\operatorname{tr}
\left(
    \left(
      \boldsymbol{I}
      -\boldsymbol{V}^*_{\rm and}(\boldsymbol{V}^\dagger_{\rm and})^*
    \right)
    \boldsymbol{X}^*
    \boldsymbol{X}
    \left(
        \boldsymbol{I}
        -\boldsymbol{V}^\dagger_{\rm and}\boldsymbol{V}_{\rm and}
    \right)
\right)
\tag{20}
\\
&=\operatorname{tr}
\left(
    \left(
        \boldsymbol{I}
        -\boldsymbol{P}^*
    \right)
    \boldsymbol{S}
    \left(
        \boldsymbol{I}
        -\boldsymbol{P}
    \right)
\right)
\tag{21}
\\
&=\operatorname{tr}(\boldsymbol{S})
-\operatorname{tr}(\boldsymbol{SP})
-\operatorname{tr}(\boldsymbol{P}^*\boldsymbol{S})
-\operatorname{tr}(\boldsymbol{P}^*\boldsymbol{SP})
\tag{22}
\\
&=
\operatorname{tr}(\boldsymbol{S})
-2\operatorname{tr}(\boldsymbol{SP})
-\operatorname{tr}(\boldsymbol{P}^*\boldsymbol{SP})
\quad
\because\boldsymbol{S}=\boldsymbol{S}^*
\tag{23}
\end{align}

この最適化問題は解析的に求めることが出来ないため、準ニュートン法(BFGS)のような数値最適化アルゴリズムを採用する。

実装例

Optimized DMDは数値最適化によって最小化問題を解くので、その初期値としてExact DMDを用いる。
まずは、元のシステム再構築のためにVandermonde行列を生成する関数を作成：

# Vandermonde matrix
def vandermonde(mu, ts):
    mu = mu.reshape(-1, 1)
    ts = ts.reshape(1, -1)
    return np.exp(mu @ ts)

これを用いてExact DMDを行う関数を作成：

# Exact DMD
def dmd(X, rank, dt=1.0):
    # Split into _X and Xprime
    Xprime, _X = X[1:, 1:], X[:, :-1]

    # 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
    eigs, W = np.linalg.eig(Atilde)
    mus = np.log(eigs) / dt

    # Reconstruction
    Phi = Xprime @ np.linalg.solve(Sr.T, Vhr).T @ W @ np.diag(1 / eigs)

    x1tilde = Sr @ Vhr[:, 0]
    b = np.linalg.solve(W @ np.diag(eigs), x1tilde)

    # Setup time step and Vandermonde matrix
    ts = np.arange(X.shape[1]) * dt
    V = vandermonde(mus, ts)

    # Reconstruct the DMD modes
    Xr = Phi @ np.diag(b) @ V

    return (Phi, eigs, mus, b, Xr)

複素数の最適化問題となるので、複素数セット $\lbrace\mu_1,\mu_2,\cdots,\mu_r\rbrace$ を実数セット $\lbrace a_1,a_2,\cdots,\ a_r,b_1,b_2,\cdots,b_r\rbrace$ のように変形してから最適化アルゴリズムをかける。
よって、複素数セット↔︎実数セットの変換関数を作成する：

# Convert the complex array to a real array
def mu2ab(mu):
    mu = np.vstack([np.real(mu), np.imag(mu)])
    ab = mu.flatten()
    return ab

# Convert the real array back to a complex array
def ab2mu(ab, N):
    r = len(ab) // (2 * N)
    ab = ab.reshape(-1, N)
    mu = ab[:r] + 1j * ab[r:]
    return mu

続いて、式 $(22)$ を再現した損失とその勾配（ヤコビアン）を計算する関数を作成する。勾配の導出ステップについては次の通り。$\boldsymbol{\mu} := \lbrace\mu_1,\mu_2,\cdots,\mu_r\rbrace$ と置いて、

損失 $f$ を行列 $\boldsymbol{V}_{\rm and}$ で微分
行列 $\boldsymbol{V}_{\rm and}$ を連続固有値 $\boldsymbol{\mu}$ で微分
チェインルールで勾配 $\nabla f$ を求める

1. 損失 $f$ を行列 $\boldsymbol{V}_{\rm and}$ で微分

\begin{align}
f(\boldsymbol{\mu})
&=
\operatorname{tr}(\boldsymbol{S})
-2\operatorname{tr}(\boldsymbol{SP})
-\operatorname{tr}(\boldsymbol{P}^*\boldsymbol{SP})
\tag{24}
\\
\Rightarrow
\cfrac{
    \partial f
}{
    \partial\boldsymbol{V}_{\rm and}
}
&=
-2\,
\cfrac{
    \partial
}{
    \partial\boldsymbol{V}_{\rm and}
}
\operatorname{tr}(\boldsymbol{SP})
+
\cfrac{
    \partial
}{
    \partial\boldsymbol{V}_{\rm and}
}
\operatorname{tr}
(
    \boldsymbol{P}^*
    \boldsymbol{SP}
)
\tag{25}
\end{align}

ここで各項の偏微分について考えるが、$\boldsymbol{V}_{\rm and}^\dagger$ を定数と見做し、$\boldsymbol{V}_{\rm and}^\dagger\boldsymbol{V}_{\rm and}$ の微小変化は $\boldsymbol{V}_{\rm and}$ のみによって支配されるものと考える。

$\operatorname{tr}(\boldsymbol{SP})$の微分

\cfrac{
    \partial
}{
    \partial\boldsymbol{V}_{\rm and}
}
\operatorname{tr}(\boldsymbol{SP})
=
\cfrac{
    \partial
}{
    \partial\boldsymbol{V}_{\rm and}
}
\operatorname{tr}(
    \boldsymbol{SV}^\dagger_{\rm and}
    \boldsymbol{V}_{\rm and}
)
=
\boldsymbol{V}^{\dagger *}_{\rm and}
\boldsymbol{S}
\tag{26}

$\operatorname{tr}(\boldsymbol{P}^*\boldsymbol{SP})$の微分

\begin{align}
\cfrac{
    \partial
}{
    \partial\boldsymbol{V}_{\rm and}
}
\operatorname{tr}
(
    \boldsymbol{P}^*
    \boldsymbol{SP}
)
&=
\cfrac{
    \partial
}{
    \partial\boldsymbol{V}_{\rm and}
}
\operatorname{tr}
(
    \boldsymbol{V}^*_{\rm and}
    \boldsymbol{V}^{\dagger *}_{\rm and}
    \boldsymbol{SV}^\dagger_{\rm and}
    \boldsymbol{V}_{\rm and}
)
\tag{27}
\\
&=
\left(
    \boldsymbol{V}^{\dagger *}_{\rm and}
    \boldsymbol{SV}^\dagger_{\rm and}
\right)
\boldsymbol{V}_{\rm and}
+
\left(
   \boldsymbol{V}^{\dagger *}_{\rm and}
    \boldsymbol{SV}^\dagger_{\rm and}
\right)^*
\boldsymbol{V}_{\rm and}
\tag{28}
\\
&=
2
\left(
    \boldsymbol{V}^{\dagger *}_{\rm and}
    \boldsymbol{SV}^\dagger_{\rm and}
\right)
\boldsymbol{V}_{\rm and}
\tag{29}
\end{align}

故に式 $(25)$ は：

\cfrac{
    \partial f
}{
    \partial\boldsymbol{V}_{\rm and}
}
=
2
\boldsymbol{V}^{\dagger *}_{\rm and}
\boldsymbol{S}
\left(
    \boldsymbol{P}-\boldsymbol{I}
\right)
\tag{30}

2. 行列 $\boldsymbol{V}_{\rm and}$ を連続固有値 $\boldsymbol{\mu}$ で微分

式 $(10)$ のVandermonde行列は離散時間となっているが、実際には $\mu_j:=\log{\lambda_j}/\Delta t$ とし、連続時間を考える：

\begin{align}
\boldsymbol{V}_{\rm and}
&=
\begin{bmatrix}
    e^{\mu_1t_1}  & \cdots & e^{\mu_1t_\tau} \\
    \vdots & \ddots & \vdots \\
    e^{\mu_rt_1} & \cdots & e^{\mu_rt_\tau}
\end{bmatrix}
\tag{31}
\\
\Rightarrow
\frac{
    \partial\boldsymbol{V}_{\rm and}
}{
    \partial\mu_j
}
&=
\begin{bmatrix}
    0 & \cdots & 0 \\
    \vdots & \ddots & \vdots \\
    0 & \cdots & 0 \\
    t_1e^{\mu_jt_1} & \cdots & t_me^{\mu_jt_\tau} \\
    0 & \cdots & 0 \\
    \vdots & \ddots & \vdots \\
    0 & \cdots & 0 \\
\end{bmatrix}
\tag{32}
\end{align}

$\frac{\partial\boldsymbol{V}_{\rm and}}{\partial\boldsymbol{\mu}}$ はサイズ $\mathbb{R}^{r\times\tau\times r}$ のテンソルとなり扱いが煩雑なので、縮約記法を用いる：

\begin{align}
\cfrac{
    \partial\boldsymbol{V}_{{\rm and}\,i,j}
}{
    \partial\mu_i
}
&=
t_je^{\mu_it_j}
\tag{33}
\\
\Rightarrow
\left(
   \cfrac{
       \partial\boldsymbol{V}_{\rm and}
   }{
       \partial\mu_i
   }
\right)_j
&=
\left(
   \boldsymbol{T}\odot\boldsymbol{V}_{\rm and}
\right)_{i,j}
\tag{34}
\end{align}

3. チェインルールで勾配 $\nabla f$ を求める

今得られた2つの導関数を組み合わせ、合成関数のチェインルールを用いて最終的に欲しい勾配を計算する。
Wirtinger 微分に基づき、実スカラー関数 $f$ の勾配は $\nabla f$は：

\begin{align}
\nabla f
=
2\cdot
\cfrac{
    \partial f
}{
    \partial
    \overline{\mu_i}
}
&=
2\cdot
\overline{
   \left(
      \cfrac{
          \partial f
      }{
          \partial
          \mu_i
      }
   \right)
}
\tag{35}
\\
&=
\sum_{j=1}^\tau
\overline{
   \left(
      \frac{
          \partial f
      }{
          \partial\boldsymbol{V}_{{\rm and}\,i,j}
      }
   \right)
}
\cdot
\overline{
   \left(
      \cfrac{
          \partial\boldsymbol{V}_{{\rm and}\,i,j}
      }{
          \partial\mu_i
      }
   \right)
}
\tag{36}
\\
&=
\sum_{j=1}^\tau
\overline{
   \left(
      \frac{
          \partial f
      }{
          \partial\boldsymbol{V}_{\rm and}
      }
   \right)_{i,j}
}
\cdot
\overline{
   t_je^{\mu_it_j}
}
\tag{37}
\end{align}

これを実装したクラスが下記となる：

class DMDLoss:
    def __init__(self, S: np.ndarray, ts: np.ndarray, gamma: float = 1e-3):
        self.S = S
        self.ts = ts
        self.gamma = gamma
        self._cache = {}

    def evaluate(self, ab: np.ndarray) -> tuple[float, np.ndarray]:
        key = ab.tobytes()
        if key in self._cache:
            return self._cache[key]

        mu = ab2mu(ab, 1)
        V = vandermonde(mu, self.ts)
        pinvV = np.linalg.pinv(V)
        P = pinvV @ V

        S = self.S
        f = np.real(np.trace(S - 2 * S @ P + P.conj().T @ S @ P))

        df_dV = 2 * pinvV.conj().T @ S @ (P - np.eye(len(self.ts)))
        dV_dmu = V * self.ts
        df_dmu = np.sum(df_dV.conj() * dV_dmu.conj(), axis=1)
        g = 2 * mu2ab(df_dmu)

        # L2 regularization
        reg = self.gamma / len(ab) / 2
        f += reg * np.dot(ab[: len(ab) // 2], ab[: len(ab) // 2])
        g[: len(g) // 2] += reg * 2 * ab[: len(ab) // 2]

        self._cache[key] = (f, g)
        return f, g

    def fun(self, ab: np.ndarray) -> float:
        return self.evaluate(ab)[0]

    def grad(self, ab: np.ndarray) -> np.ndarray:
        return self.evaluate(ab)[1]

最後にモード振幅 $\boldsymbol{b}$ について、Exact DMDでは以下の様に与えられていた：

\begin{align}
\begin{cases}
\tilde{\boldsymbol{x}}_1
=
\left(
   \tilde{\boldsymbol{\Sigma}}
   \tilde{\boldsymbol{V}}
\right)_{:,1}\\
\boldsymbol{b}
=
\left(
    \boldsymbol{W}
    \boldsymbol{\Lambda}
\right)^{-1}
\tilde{\boldsymbol{x}}_1
\end{cases}
\tag{38}
\end{align}

行列 $\boldsymbol{W}$ と行列 $\boldsymbol{\Lambda}$ はそれぞれ線形作用素 $\boldsymbol{A}$ の固有行列と固有値行列である。この定義に倣えば、Optimized DMDにおけるモード振幅 $\boldsymbol{b}$ は次の様になる：

\boldsymbol{x}_1
\simeq\sum_{j=1}^r
\psi_j
\tag{39}

各モードに対応する振幅を $b_j:=1/|\psi_j|_2$ とすることで：

\boldsymbol{x}_1
\simeq\sum_{j=1}^r
\psi_j
=
\sum_{j=1}^r
\cfrac{
    1
}{
    \|\psi_j\|_2
}
\cdot
\hat{\psi}_j
\tag{40}

故に、Optimized DMDにおけるモード振幅は $b_j:=1/|\psi_j|_2=\operatorname{tr}(\boldsymbol{\Psi}^*\boldsymbol{\Psi})$ と定義できる。

Optimized DMD全体の関数は以下の通り：

def opt_dmd(X, rank, dt=1.0):
    # Exact DMD for mode initialization
    _, _, mus_init, _, _ = dmd(X, rank)

    # Optimization problem setup
    init_ab = mu2ab(mus_init).reshape(2 * rank)

    # Setup time step
    ts = np.arange(X.shape[1]) * dt

    # Obtain loss function and gradient function of it
    S = X.T @ X
    loss = DMDLoss(S, ts)

    # Boundary conditions for real parts
    bounds = [
        (-1e10, 0.1) if i < len(init_ab) // 2 else (None, None)
        for i in range(len(init_ab))
    ]

    # Solve optimization problem
    res = minimize(
        fun=loss.fun,
        jac=loss.grad,
        x0=init_ab,
        bounds=bounds,
        method="L-BFGS-B",
        options={"maxiter": 1000, "ftol": 1e-9, "disp": True},
    )
    mus = ab2mu(res.x, 1)
    eigs = np.exp(mus * dt).flatten()

    # DMD mode calculation
    V = vandermonde(mus, ts)
    W = X @ np.linalg.pinv(V)

    # Calculation of mode amplitude
    b = np.diag(W.conj().T @ W)

    # Reconstruction
    Xr = W @ V
    return (W, eigs, mus, b, Xr)

Optimized DMDの実装は以上となるので、Exact DMDと再構築誤差比較を行ってみる。時間ステップごとに状態の誤差を算出する関数を作成：

# Calculate the reconstruction error percentage
def reconstruction_error(func, *args, **kwargs):
    _, _, _, _, Xr = func(*args, **kwargs)
    return np.linalg.norm(X - Xr.real, axis=0) * 100 / np.linalg.norm(X, axis=0)

これを用いて誤差をプロット：

# Measure the accuracy of the exact DMD and optimized DMD
exact_errors = reconstruction_error(dmd, X, 20)
opt_errors = reconstruction_error(opt_dmd, X, 20)

# Plot the errors
plt.figure(figsize=(10, 6))
plt.plot(np.arange(X.shape[1]), exact_errors, label="Exact DMD", marker="s", linestyle="--")
plt.plot(np.arange(X.shape[1]), opt_errors, label="Optimized DMD", marker="o", linestyle="-")

# Customize the plot
plt.xlabel("Time Step")
plt.ylabel("Error (%)")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.6)

# Display the plot
plt.show()

[出力結果]

Exact DMDと比べて復元誤差がかなり改善された。
因みに、Optimized DMDだけの誤差を描画すると下記の様になる。

Sparsity-promoting DMD (M.R.Jovanovic et al., 2013)

Exact DMDでは得られる全てのモードを用いてデータを再構成するが、Sparsity-promoting DMD（以下、SpDMD）は少数の支配的なDMDモードで元データを再現することで、モデル解釈性や計算効率向上を図る手法である。

手法と数式モデル

SpDMDの基本的な考え方は次の通り。まず、Exact DMDと同様に連続時間系のシステム再構築を考える：

\boldsymbol{X}
\simeq
\boldsymbol{\Phi}{\rm diag}
\left(
    \boldsymbol{b}
\right)
\boldsymbol{V}_{\rm and}
\tag{41}

Exact DMDではモード振幅 $\boldsymbol{b}$ を最初のスナップショットの縮約表現 $\tilde{\boldsymbol{x}}_1$ を用いて推定する：

\boldsymbol{b}
=
\left(
    \boldsymbol{W}
    \boldsymbol{\Lambda}
\right)^{-1}
\tilde{\boldsymbol{x}}_1
\tag{42}

行列 $\boldsymbol{W}$ と行列 $\boldsymbol{\Lambda}$ はそれぞれ線形作用素 $\boldsymbol{A}$ の固有行列と固有値行列である。
この問題点は、モード推定が初期状態 $\boldsymbol{x}_1$ に依存しており、システム全体に対して最適化されていない事である。前章のOptimized DMDでは固有値 $\lambda_j$ を推定していたが、SpDMDではモード振幅 $\boldsymbol{b}$ のみを最適化することを図る：

\min_{b_j}
\|
    \boldsymbol{X}
    -
    \boldsymbol{\Phi}
    {\rm diag}
   \left(
       \boldsymbol{b}
   \right)
   \boldsymbol{V}_{\rm and}
\|_F^2
\tag{43}

さらに、ここに対して$\ell_1$正規化項を追加して $b_j$ のスパース性を促す：

\begin{align}
J(\boldsymbol{b})
:=&
\|
    \boldsymbol{X}
    -
    \boldsymbol{\Phi}
    {\rm diag}
   \left(
       \boldsymbol{b}
   \right)
   \boldsymbol{V}_{\rm and}
\|_F^2
+
\gamma
\|
\boldsymbol{b}
\|_1
\tag{44}
\\
&\Rightarrow
\min_{\boldsymbol{b}}
J(\boldsymbol{b})
\tag{45}
\end{align}

$\gamma$ はスパース性制御パラメータであり、これが大きいほど少数のモードが選択される。
さて、ここに縮約表現 $\boldsymbol{X}\simeq\tilde{\boldsymbol{U}}\tilde{\boldsymbol{\Sigma}}^{-1}\tilde{\boldsymbol{V}}^*$ 及び、射影モードの定義 $\boldsymbol{\Phi}\simeq\tilde{\boldsymbol{U}}\boldsymbol{W}$ を代入すると：

\begin{align}
J(\boldsymbol{b})
&=
\|
    \tilde{\boldsymbol{U}}
    \tilde{\boldsymbol{\Sigma}}^{-1}
    \tilde{\boldsymbol{V}}^*
    -
    \tilde{\boldsymbol{U}}
    \boldsymbol{W}
    {\rm diag}
    \left(
        \boldsymbol{b}
    \right)
    \boldsymbol{V}_{\rm and}
\|_F^2
+
\gamma
\|
\boldsymbol{b}
\|_1
\tag{46}
\\
&
\equiv
\|
    \tilde{\boldsymbol{\Sigma}}^{-1}
    \tilde{\boldsymbol{V}}^*
    -
    \boldsymbol{W}
    {\rm diag}
    \left(
        \boldsymbol{b}
    \right)
    \boldsymbol{V}_{\rm and}
\|_F^2
+
\gamma
\|
\boldsymbol{b}
\|_1
\tag{47}
\end{align}

実は、式 $(47)$ の第一項目は $\boldsymbol{b}$ で微分可能であり、凸最適化であることから第二項がなければ、解析的に最適解を求めることができる。

実際に、第一項目の実スカラー関数 $\tilde{J}(\boldsymbol{b})=|
\tilde{\boldsymbol{\Sigma}}^{-1}
\tilde{\boldsymbol{V}}^*
-
\boldsymbol{W}
{\rm diag}
\left(
\boldsymbol{b}
\right)
\boldsymbol{V}_{\rm and}
|_F^2\in\mathbb{R}$ の最小値を得るには、Wirtinger微分における $\partial\tilde{J}/\partial\overline{\boldsymbol{b}}=0$ を解けば良い。$\boldsymbol{Y}
:=\tilde{\boldsymbol{\Sigma}}^{-1}\tilde{\boldsymbol{V}}^*
=\tilde{\boldsymbol{U}}^*\boldsymbol{X}$、$\boldsymbol{B}:={\rm diag}\left(\boldsymbol{b}\right)$ と置いて：

\begin{align}
&&
\tilde{J}
(
    \boldsymbol{b}
)
&=
\|
    \boldsymbol{Y}
    -
    \boldsymbol{W}
    \boldsymbol{B}
    \boldsymbol{V}_{\rm and}
\|_F^2
\tag{48}
\\
&&
&=
\operatorname{tr}
\left[
    \left(
        \boldsymbol{Y}
        -
        \boldsymbol{W}  
        \boldsymbol{B}
        \boldsymbol{V}_{\rm and}
    \right)
     \left(
        \boldsymbol{Y}
        -
        \boldsymbol{W}
        \boldsymbol{B}
        \boldsymbol{V}_{\rm and}
    \right)^*
\right]
\tag{49}
\\
&&
&=
\operatorname{tr}
\left(
    \boldsymbol{Y}
    \boldsymbol{Y}^*
\right)
-
\operatorname{tr}
\left(
    \boldsymbol{W}
    \boldsymbol{B}
    \boldsymbol{V}_{\rm and}
    \boldsymbol{Y}^*
\right)
-
\operatorname{tr}
\left(
    \boldsymbol{Y}
    \boldsymbol{V}_{\rm and}^*
    \boldsymbol{B}^*
    \boldsymbol{W}^*
\right)
+
\operatorname{tr}
\left(
    \boldsymbol{W}
    \boldsymbol{B}
    \boldsymbol{V}_{\rm and}
    \boldsymbol{V}_{\rm and}^*
    \boldsymbol{B}^*
    \boldsymbol{W}^*
\right)
\tag{50}
\\
&\Rightarrow&
\cfrac{
    \partial\tilde{J}
}{
    \partial
    \boldsymbol{B}^*
}
&=
-
\boldsymbol{V}_{\rm and}
\boldsymbol{Y}^*
\boldsymbol{W}
+
\boldsymbol{V}_{\rm and}
\boldsymbol{V}_{\rm and}^*
\boldsymbol{B}^*
\boldsymbol{W}^*
\boldsymbol{W}
=0
\tag{51}
\\
&\Rightarrow&
&
\boldsymbol{W}^*
\boldsymbol{Y}
\boldsymbol{V}_{\rm and}^*
=
\boldsymbol{W}^*
\boldsymbol{W}
\boldsymbol{B}
\boldsymbol{V}_{\rm and}
\boldsymbol{V}_{\rm and}^*
\tag{52}
\end{align}

ここで、$\boldsymbol{B}$ が対角行列であることから対角成分 $j=l$ に注目して：

\begin{align}
&&
\left(
    \boldsymbol{W}^*
    \boldsymbol{Y}
    \boldsymbol{V}_{\rm and}^*
\right)_{j,j}
=
\sum_k
b_k
\left(
    \boldsymbol{W}^*
    \boldsymbol{W}
\right)_{j,k}
\left(
    \boldsymbol{V}_{\rm and}
    \boldsymbol{V}_{\rm and}^*
\right)_{k,j}
\tag{53}
\\
&\Rightarrow&
\text{diag}
\left(
    \boldsymbol{W}^*
    \boldsymbol{Y}
    \boldsymbol{V}_{\rm and}^*
\right)
=
\left(
    \boldsymbol{W}^*
    \boldsymbol{W}
\right)
\odot
\left(
    \boldsymbol{V}_{\rm and}
    \boldsymbol{V}_{\rm and}^*
\right)^\top
\boldsymbol{b}
\tag{54}
\\
&\Rightarrow&
\boldsymbol{b}
=
\left(
    \left(
        \boldsymbol{W}^*
        \boldsymbol{W}
    \right)
    \odot
    \left(
        \boldsymbol{V}_{\rm and}
        \boldsymbol{V}_{\rm and}^*
    \right)^\top
\right)^{-1}
\text{diag}
\left(
    \boldsymbol{W}^*
    \boldsymbol{Y}
    \boldsymbol{V}_{\rm and}^*
\right)
\tag{55}
\end{align}

本手法では $\ell_1$ 正則化項がついているため、ADMM（Alternating Direction Method of Multipliers）を用いて最適解を近似する。

ADMMによる求解

元の最適化問題を変数分離して処理をする：

\min_{\boldsymbol{b},\boldsymbol{\alpha}}
J(\boldsymbol{b})
+
\gamma
\|
\alpha
\|_1
\quad
\text{subject to}
\quad
\boldsymbol{b}-\boldsymbol{\alpha}
=
\boldsymbol{0}
\tag{56}

これは式 $(47)$ と等価である。
ここで、ラグランジュ関数法を用いる。ラグランジュ乗数を $\boldsymbol{\zeta}=\lbrace\zeta_j\rbrace$ と置いて：

\begin{align}
\mathcal{L}
(
    \boldsymbol{b},
    \boldsymbol{\alpha},
    \boldsymbol{\zeta}
)
&=
J(\boldsymbol{b})
+\gamma
\|
    \alpha
\|_1
+
\sum_{j=1}^r
\Re{
    \left(
        \zeta_j(
            b_j-\alpha_j
        )
    \right)
}
\tag{57}
\\
&=
J(\boldsymbol{b})
+\gamma
\|
    \alpha
\|_1
+
\Re{
    \left(
        \boldsymbol{\zeta}^*
        (
            \boldsymbol{b}
            -\boldsymbol{\alpha}
        )
    \right)
}
\tag{58}
\end{align}

しかしながら、このままでは数値安定性や収束性が悪いので、ペナルティ項 $\frac{1}{2}\rho|\boldsymbol{b}-\boldsymbol{\alpha}|_2^2$ を加え、拡張ラグランジュ関数法に切り替える：

\begin{align}
\mathcal{L}_\rho
(
    \boldsymbol{b},
    \boldsymbol{\alpha},
    \boldsymbol{\zeta}
)
&=
J(\boldsymbol{b})
+\gamma
\|
    \alpha
\|_1
+
\Re{
    \left(
        \boldsymbol{\zeta}^*
        (
            \boldsymbol{b}
            -\boldsymbol{\alpha}
        )
    \right)
}
+
\frac{1}{2}
\rho
\|
    \boldsymbol{b}-\boldsymbol{\alpha}
\|_2^2
\tag{59}
\\
&=
J(\boldsymbol{b})
+\gamma
\|
    \alpha
\|_1
+
\cfrac{1}{2}
\lbrace
    \boldsymbol{\zeta}^*
    (
        \boldsymbol{b}
        -\boldsymbol{\alpha}
    )
    +
    (
        \boldsymbol{b}^*
        -\boldsymbol{\alpha}^*
    )
    \boldsymbol{\zeta}
\rbrace
+
\frac{1}{2}
\rho
\|
    \boldsymbol{b}-\boldsymbol{\alpha}
\|_2^2
\tag{60}
\\
&=
J(\boldsymbol{b})
+\gamma
\|
    \alpha
\|_1
+
\cfrac{1}{2}
\lbrace
    \boldsymbol{\zeta}^*
    (
        \boldsymbol{b}
        -\boldsymbol{\alpha}
    )
    +
    (
        \boldsymbol{b}^*
        -\boldsymbol{\alpha}^*
    )
    \boldsymbol{\zeta}
   +
   \rho
   \|
       \boldsymbol{b}-\boldsymbol{\alpha}
   \|_2^2
\rbrace
\tag{61}
\end{align}

各反復 $h$ において $\rho$ を更新し（更新しなくても収束はする）、$\boldsymbol{\zeta}_j$ も以下の式によって更新される:

\boldsymbol{\zeta}_j^{(h+1)}
\leftarrow
\boldsymbol{\zeta}_j^{(h)}
+
\rho
\left(
    \boldsymbol{b}^{(h+1)}
    -\boldsymbol{\alpha}^{(h+1)}
\right)
\tag{62}

ただし、$\boldsymbol{b}^{(h+1)}$ は $h$ 番目の反復における問題を解いた解 $\text{arg min}_{\boldsymbol{b}}{\mathcal{L}_\rho^{(h)}}$ である。$\boldsymbol{\alpha}^{(h+1)}$ についても同じ。
停止条件は以下の2つを設ける：

\begin{align}
\text{制約違反：}&
\|
    \boldsymbol{b}^{(h+1)}
    -\boldsymbol{\alpha}^{(h+1)}
\|_2
&\leq
\varepsilon_{\rm prim}
\tag{63}
\\
\text{制約差分：}&
\|
    \boldsymbol{\alpha}^{(h+1)}
    -\boldsymbol{\alpha}^{(h)}
\|_2
&\leq
\varepsilon_{\rm dual}
\tag{64}
\end{align}

実装例

このspDMDの実装はとても遅く、同僚から違うのではないかと指摘されていますので、あまり参考にしないでください。アルゴリズムの解釈なのか実装が原因なのかもまだ調べられていません。

はじめに、式 $(55)$ の解析ソルバーを実装：

# amplitude spectrum solver function
def solve_b(Y, W, V, lam=0):
    # Inner products
    WW = W.conj().T @ W
    VV = V @ V.conj().T

    # Hadamard product
    P = WW * VV.T
    cond_P = np.linalg.cond(P)
    print(f"Condition number of P: {cond_P:.2e}")

    # q vector: diag(Vand A* W)
    Q = W.conj().T @ Y @ V.conj().T
    q = np.diag(Q)
    
    # Solve
    if lam:
        P = P + lam * np.linalg.norm(P, 2) * np.eye(P.shape[0])

    return np.linalg.solve(P, q)

もし、行列 $\boldsymbol{P}$ の条件数が大きい場合には、変数 lam で正則化項を入れて調整すると良い。
そして、これを使ったモード振幅最適化DMD b_opt_DMD：

# DMD w/ the analytical mode amplitude
def b_opt_dmd(X, rank, dt=1.0):
    # Split into _X and Xprime
    Xprime, _X = X[:, 1:], X[:, :-1]

    # 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
    eigs, W = np.linalg.eig(Atilde)
    mus = np.log(eigs) / dt

    # Reconstruction
    Phi = Xprime @ np.linalg.solve(Sr.T, Vhr).T @ W @ np.diag(1 / eigs)

    # Setup time step and Vandermonde matrix
    ts = np.arange(X.shape[1]) * dt
    V = vandermonde(mus, ts)

    # Calculate the mode amplitude
    Y = Ur.conj().T @ X
    b = solve_b(Y, W, V)

    # Reconstruct the DMD modes
    Xr = Phi @ np.diag(b) @ V

    return (Phi, eigs, mus, b, Xr)

さて、以下から $\ell_1$ 正則化項付きSpDMDの実装に入る。まずは、目的関数 $J(\boldsymbol{b})$ ：

def obj(ab, X, Phi, V, gamma):
    b = ab2mu(ab, 1).squeeze()
    return np.linalg.norm(
        X - Phi @ np.diag(b) @ V, ord="fro"
    ) ** 2 + gamma * np.linalg.norm(b, ord=1)

続いて、拡張ラグランジュ関数 $\mathcal{L}_\rho(\boldsymbol{b},\boldsymbol{\alpha},\boldsymbol{\zeta})$：

def aug_lagrangian(ab, alpha, X, Phi, V, zeta, rho, gamma):
    diff = ab2mu(ab - alpha, 1)
    return (
        obj(ab, X, Phi, V, gamma)
        +np.real(
            ab2mu(zeta, 1).conj().T @ diff
        )
        + 1/2 * (
            + rho * np.linalg.norm(diff, ord=2) ** 2
        )
    )

ソフト閾値処理。これによってスパース化を実現する：

def soft_thresholding(x, thresh):
    return np.maximum(0, 1 - thresh / np.abs(x)) * x

SpDMD全体の関数は以下の通り：

def sp_dmd(X, rank, dt=1.0, rho=1.0, gamma=4e3, tol=1e-3, max_iter=10):
    # Exact DMD for mode initialization
    Phi, eigs, mus, b_init, _ = dmd(X, rank)
    
    # Optimization problem setup
    b = mu2ab(b_init)
    alpha = np.copy(b)
    zeta = np.zeros_like(b)

    # Setup Vandermonde matrix
    ts = np.arange(X.shape[1]) * dt
    V = vandermonde(mus, ts)

    # AMDD algorithm
    for _ in tqdm(range(int(max_iter))):
        fun = lambda x: aug_lagrangian(x, alpha, X, Phi, V, zeta, rho, gamma)
        res = minimize(
            fun=fun,
            x0=b,
            method="L-BFGS-B",
            options={"ftol": 1e-6, "maxiter": 100},
        )
        b = res.x
        
        # alpha-update via soft thresholding
        p_alpha = alpha.copy()
        alpha = soft_thresholding(b + zeta / rho, gamma / rho)
        
        # zeta-update (dual ascent)
        zeta += rho * (b - alpha)

        # Convergence check
        r_primal = np.linalg.norm(b - alpha)
        r_dual = np.linalg.norm(alpha - p_alpha)
        if r_primal < tol and r_dual < tol:
            break
        
    # Reconstruct the DMD modes
    b = ab2mu(res.x, 1).squeeze()
    Xr = Phi @ np.diag(b) @ V
    
    return (Phi, eigs, mus, b, Xr)

ただし、スパース性制御パラメータ $\gamma$ は、データに依存するので、 $\tilde{J}$ と $|\boldsymbol{b}|_1$ の大きさを比較して決定した。
先ほど同様に、Exact DMDと並べて時間ごと復元誤差をプロットすると以下のようになる。

また、b-opt-DMDとSpDMDそれぞれのモード振幅と周波数を標準出力すると次の通りなる：

======================================================================
Index 絶対固有値       周波数 (Hz)    b-opt-DMDモード振幅   SpDMDモード振幅      
----------------------------------------------------------------------
1     0.484637       0.500000       3.416869            0.000002       
2     1.000002       0.297733       0.565270            0.000001       
3     1.000002       -0.297733      0.565270            0.000001       
4     1.000030       0.264640       1.010828            0.000001       
5     1.000030       -0.264640      1.010828            0.000001       
6     0.999991       0.231557       1.781063            0.000003       
7     0.999991       -0.231557      1.781063            0.000004       
8     0.999994       0.198477       3.050829            0.000001       
9     0.999994       -0.198477      3.050829            0.000009       
10    1.000002       0.165397       6.274135            0.000017       
11    1.000002       -0.165397      6.274135            0.000008       
12    1.000001       0.132317       11.849283           0.110507       
13    1.000001       -0.132317      11.849283           0.110017       
14    1.000000       0.000000       290.832289          278.770667     
15    1.000000       0.099238       29.458344           17.226386      
16    1.000000       -0.099238      29.458344           17.229267      
17    1.000000       0.033079       104.819063          92.366314      
18    1.000000       -0.033079      104.819063          92.365891      
19    1.000000       0.066158       37.849752           25.533114      
20    1.000000       -0.066158      37.849752           25.533130      
======================================================================

モード 14 は平均場なので、モード 15 ~ 20 が優位なモードであると示唆している。

Randomized DMD (N.B.Erichson et al., 2017)

流体などの実データだと、スナップショット行列 $\boldsymbol{X}\in\mathbb{R}^{m\times\tau}$ は $m\gg\tau$ であり、メモリなど必要となる計算リソースが膨大になる。Randomized DMD（以下、rDMD）は、ランダム射影で列空間を近似し、低次元空間でのDMDを実現する手法である。

rDMDのイメージ図。高次元のスナップショットデータ $\boldsymbol{X},,\boldsymbol{X}'$ から低次元のスナップショット行列 $\boldsymbol{B},,\boldsymbol{B}'$ を求める。これらを用いて近似DMDモード $\tilde{\boldsymbol{\Phi}}$ と固有値 $\boldsymbol{\Lambda}$ を計算し、元のDMDモードを復元する。

手法と数式モデル

rDMDは、1.次元を乱択的に圧縮、2.DMDを実行、の2段階の確率的手法に基づいている。この確率的フレームワーク（Probabilistic Framework）はHalko et al. が組み立てたものであり、以前に紹介したRandomized SVDにも用いられている。rDMDはランダム低ランク近似アルゴリズムをDMDに移植したものである。

近似列空間基底 $\boldsymbol{Q}$ の構築

はじめに、入力となる行列 $\boldsymbol{X}\in\mathbb{R}^{m\times\tau}$ の列空間を低次元で近似するため、その列空間を張る直交基底の近似行列 $\boldsymbol{Q}$ を構築する。

ランダム射影

階数 $r\leq\tau$ を決め、入力スナップショット行列 $\boldsymbol{X}\in\mathbb{R}^{m\times\tau}$ に対し、列数 $l=r+p, p\in\mathbb{N}$ の乱数行列 $\boldsymbol{\Omega}\in\mathbb{R}^{\tau\times l}$ を掛ける：

\boldsymbol{Y}
=
\boldsymbol{X\Omega}
\tag{65}

通常、乱数行列 $\boldsymbol{\Omega}$ は独立ガウス分布から生成される。
このようにして得られた行列 $\boldsymbol{Y}$ の列空間は高確率で $\boldsymbol{X}$ の支配的 $r$ 次元部分空間を捉える。

直交化

$\boldsymbol{Y}$ に対して低ランクQR分解 $\boldsymbol{Y}=\boldsymbol{QR}$ を行い、列直交基底 $\boldsymbol{Q}\in\mathbb{R}^{m\times l}$ を得る：

\boldsymbol{X}
\simeq
\boldsymbol{QQ}^*
\boldsymbol{X}
\tag{66}

オーバーサンプリング

実は、乱数行列 $\boldsymbol{\Omega}$ の列数を $l=r$ ではなく、$p$ だけ（通常5~10程度）増やすには理由がある。
ランダム射影を用いて低ランク近似を作るとき、理想的には：

\text{rank}
\left(
    \boldsymbol{X}
\right)
=
r
\tag{67}

であり、支配的な $r$ 個の特異値 $\lbrace\sigma_i\rbrace_{i=1}^r$ だけが非ゼロで、残りの特異値 $\lbrace\sigma_i\rbrace_{i>r}$ は完全ゼロであることが望ましい。
しかし、現実データでは準低ランクの状況となっており、

残りの全ての特異値 $\lbrace\sigma_i\rbrace_{i>r}$ が完全ゼロになることは極めて稀
故に、$\boldsymbol{X}$ の列空間は「ピッタリ $r$ 次元」、ではなく、「だいたい $r$ 次元」

そこで、列数を $p$ だけ増やすことで、1.列空間が漏れる可能性を下げ、2.統計的な揺らぎが減り、特異値推定が安定する、というメリットがある。

低次元スケッチ行列 B の生成

次に、$\boldsymbol{X}\in\mathbb{R}^{m\times\tau}$ を近似列空間規定 $Q\in\mathbb{R}^{m\times l}$ で直交射影したスケッチ行列 $\boldsymbol{B}\in\mathbb{R}^{l\times m}$ を作る：

\boldsymbol{B}
=
\boldsymbol{Q}^*
\boldsymbol{X}
\tag{74}

ここで、$\boldsymbol{Q}$ が直交行列であることを用いると：

\boldsymbol{X}
\simeq
\boldsymbol{QB}
\tag{75}

と、$\boldsymbol{X}$ を近似復元できる。

パワーイテレーション q

スペクトル減衰：

\sigma_1\geq\sigma_2\geq\cdots\geq\sigma_r\tag{68}

が遅い行列に対しては $\boldsymbol{X}$ を $q$ 乗反復処理すると、減衰がより早い行列 $\boldsymbol{X}^{(q)}$ が作成できる。：

\begin{align}
\boldsymbol{X}^{(q)}
&=
\left(
    \boldsymbol{XX}^*
\right)^q
\tag{69}
\\
&=
\lbrace
    \left(
        \boldsymbol{U\Sigma V}^*
    \right)
    \left(
        \boldsymbol{U\Sigma V}^*
    \right)^*
\rbrace^q
\boldsymbol{X}
\tag{70}
\\
&=
\left(
    \boldsymbol{U\Sigma}^2
    \boldsymbol{U}^*
\right)^q
\boldsymbol{X}
\tag{71}
\\
&=
\left(
    \boldsymbol{U\Sigma}^{2q}
    \boldsymbol{U}^*
\right)
\left(
    \boldsymbol{U\Sigma V}^*
\right)
\tag{72}
\\
&=
\boldsymbol{U\Sigma}^{2q+1}
\boldsymbol{V}^*
\tag{73}
\end{align}

$q=1,2$ 程度でもかなり改善される。
ただし、I/Oの処理が増えることによる計算コスト増大に留意

Halko-Marinsson-Troppの理論的誤差境界 (N.Halko et al., 2011)

ランダム射影を用いた低ランク近似に対して、以下の平均誤差上界 $\mathcal{B}_{r,p,q}(\boldsymbol{X})$ が存在する：

\begin{align}
\mathbb{E}
\,
\|
    \boldsymbol{X}
    -\boldsymbol{QQ}^*
    \boldsymbol{X}
\|_2
&\leq
\mathcal{B}_{r,p,q}(\boldsymbol{X})
\tag{74}
\\
&=
\left[
    1
    +\sqrt{
        \cfrac{
            r
        }{
            p-1
        }
    }
    +e\sqrt{
        \cfrac{
            l
        }{
            p
        }
    }
    \cdot
    \sqrt{
        \tau
        -r
    }
\right]^{\frac{1}{2q+1}}
\sigma_{r+1}
\left(
    \boldsymbol{X}
\right)
\tag{75}
\end{align}

rDMDはこの低ランク近似を用いており、上記の理論誤差境界がそのまま適用できる。
故に、rDMDは目標階数 $r$、オーバーサンプリング $p$、パワーイテレーション $q$を調整することで、DMD誤差を明示的に制御できるという性質を持つ。

実装例

# Randomized DMD
def r_dmd(X, rank, dt=1.0, p=5, q=1, seed=42):
    if seed is not None:
        np.random.seed(seed)

    # Randomized range finder
    l = rank + p
    tau = X.shape[1]
    Omega = np.random.randn(tau, l)

    # Construct the approximated basis matrix
    Y = X @ Omega
    for _ in range(q):
        Y = X @ (X.T @ Y)
    Q, _ = np.linalg.qr(Y, mode='reduced')

    # Compression
    B = Q.T @ X
    Phi_B, eigs, mus, b, _ = dmd(B, rank)

    # Lift to high dimension
    Phi = Q @ Phi_B

    # Reconstruction
    ts = np.arange(tau) * dt
    V  = vandermonde(mus, ts)
    Xr = Phi @ (b[:, None] * V)

    return Phi, eigs, mus, b, Xr

時間別誤差をExact DMDと並べてプロットすると以下の通り。

見た目上はほぼ同じであるが、再構築誤差を計算してみると、Exact DMD： 13.549936% 、Randomized DMD： 13.549983% と小数第4位まで同じ値を取る。

Randomized DMDによるシステム再構築誤差の考察

本説では、rDMD再構築誤差の理論上界の導出を考える。

誤差分解

まず、rDMDの再構築誤差 $E_r:=|\boldsymbol{X}-\boldsymbol{Q}\tilde{\boldsymbol{\Phi}}\text{diag}(\boldsymbol{b})\boldsymbol{V}_{\rm and}|_F$ を展開する。$\hat{\boldsymbol{X}}_B:=\tilde{\boldsymbol{\Phi}}\text{diag}(\boldsymbol{b})\boldsymbol{V}_{\rm and}$ として：

\begin{align}
E_r
&=
\|
    \boldsymbol{X}
    -
    \boldsymbol{Q}
    \hat{\boldsymbol{X}}_B
\|_F
\tag{76}
\\
&=
\|
    \boldsymbol{X}
    -\boldsymbol{QQ}^*
    \boldsymbol{X}
    +
    \boldsymbol{QQ}^*
    \boldsymbol{X}
    -\boldsymbol{Q}
    \hat{\boldsymbol{X}}_B
\|_F
\tag{77}
\\
&\leq
\|
    \boldsymbol{X}
    -\boldsymbol{QQ}^*
    \boldsymbol{X}
\|_F
+
\|
    \boldsymbol{QQ}^*
    \boldsymbol{X}
    -\boldsymbol{Q}
    \hat{\boldsymbol{X}}_B
\|_F
\tag{78}
\\
\end{align}

式 $(78)$ はフロベニウスノルムの劣加法性（三角不等式）から導かれるものである。
また、この第1項目は射影誤差、第2項目は低次元（Exact）DMD再構築誤差をそれぞれ表している。
さらに、$\boldsymbol{Q}$ が直交行列であることを考えれば、$|\boldsymbol{QA}|_F=|\boldsymbol{A}|_F$ となり：

E_r
\leq
E_{\rm proj}+E_{\rm red}
,\quad
\begin{cases}
    E_{\rm proj}
    =
    \|
        \boldsymbol{X}
        -\boldsymbol{QQ}^*
        \boldsymbol{X}
    \|_F
    \\
    E_{\rm red}
    =
    \|
        \boldsymbol{QQ}^*
        \boldsymbol{X}
        -\boldsymbol{Q}
        \hat{\boldsymbol{X}}_B
    \|_F
\end{cases}
\tag{79}

射影誤差の上界

スペクトルノルム $|\cdot|_2$ とフロベニウスノルム $|\cdot|_F$ の関係から：

\|\boldsymbol{X}-\boldsymbol{QQ}^*\boldsymbol{X}\|_F
\leq
\sqrt{l}
\cdot
\|\boldsymbol{X}-\boldsymbol{QQ}^*\boldsymbol{X}\|_2
\tag{80}

\Rightarrow
E_{\rm proj}
\leq
\sqrt{l}
\cdot
\mathcal{B}_{r,p,q}
\left(
    \boldsymbol{X}
\right)
\tag{81}

低次元（Exact）DMD再構築誤差

低次元空間でのExact DMD誤差 $E_{\rm red}$ は元の空間でのExact DMD誤差 $E_{e}=|\boldsymbol{X}-\boldsymbol{\Phi}\text{diag}(\boldsymbol{b})\boldsymbol{V}_{\rm and}|_F=|\boldsymbol{X}-\hat{\boldsymbol{X}}|_F$ を用いて、次のように上界が決定される：

\begin{align}
E_{\rm red}
&=
\|
    \boldsymbol{Q}
    (
        \boldsymbol{Q}^*
        \boldsymbol{X}
        -
        \hat{\boldsymbol{X}}_B
    )
\|_F
\tag{82}
\\
&=
\|
    \boldsymbol{Q}^*
    \boldsymbol{X}
    -
    \hat{\boldsymbol{X}}_B
\|_F
\tag{83}
\\
&\leq
\min_{\text{rank}(\boldsymbol{Y})\leq r}
\|
    \boldsymbol{Q}^*
    \boldsymbol{X}
    -
    \boldsymbol{Y}
\|_F
\tag{84}
\\
&\leq
\|
    \boldsymbol{Q}^*
    \boldsymbol{X}
    -
    \boldsymbol{Q}^*
    \hat{\boldsymbol{X}}
\|_F
\tag{85}
\\
&\leq
\|
    \boldsymbol{Q}^*
    (
        \boldsymbol{X}
        -
        \hat{\boldsymbol{X}}
    )
\|_F
\tag{86}
\\
&=
\|
    \boldsymbol{X}
    -
    \hat{\boldsymbol{X}}
\|_F
=
E_e
\tag{87}
\end{align}

rDMD再構築誤差

ここまでの結果を用いると、rDMDの再構築誤差は：

E_r
\leq
\sqrt{l}
\cdot
\mathcal{B}_{r,p,q}
(\boldsymbol{X})
+E_e
\tag{88}

実際に、今回のデータにおいて $E_r$ を計算してみる。
まずは、目標階数 $r=20$ としたので、$\boldsymbol{X}$ の21番目の特異値 $\sigma_{21}$ の取得：

s = np.linalg.svd(X, compute_uv=False)
s[20]

今回のデータからは 3.589987 と出た。
故に、平均誤差上界 $\mathcal{B}_{r,p,q}(\boldsymbol{X})$ は：

\mathcal{B}_{20,5,1}(\boldsymbol{X})
=
\left[
    1
    +
    \sqrt{
        \frac{
            20
        }{
            5-1
        }
    }
    +e\sqrt{
        \frac{
            20+5
        }{
            5
        }
    }
    \cdot
    \sqrt{
        151
        -20
    }
\right]
^\frac{1}{2\cdot1+1}
\cdot
3.589987
\simeq
14.99
\tag{89}

また、Exact DMD再構築誤差は $E_e\simeq 556.42$ だったので、$E_r$ の上界は：

E_r
\leq
\sqrt{20+5}
\cdot
14.99
+
556.42
=
631.37

$|\boldsymbol{X}|_F\simeq4106.41$ なので、相対誤差の上限は $631.37/4106.41\simeq15.375%$となり、先に計算したrDMDの再構築誤差 13.549983% を抑えられていることが分かる。

実は、射影誤差 $E_{\rm proj}$ はわずか 0.51 程度であり、rDMD再構築誤差は殆ど低次元DMD再構築誤差が持っていることになる。

おわりに

本意は半年前に

Hankel DMD / Delay DMD
Extended DMD
Koopman Autoencoder (Neural DMD)

も併せて投稿したかったが、時間があまりないので執筆途中であるが本稿を投稿することにした。
上の項目も派生たち Ver.2 として調査したい。