単変量時系列データへのDMDの適用と結果の見方

# はじめに
こんにちは、(株)日立製作所[Lumada Data Science Lab.](https://www.hitachi.co.jp/products/it/lumada/about/ai/ldsl/index.html)の田中と申します。

早速ですが、電流や振動、温度、圧力などのセンサデータから異常を検知したり、異常の予兆を検知するというタスクが世の中にはあります。このとき、扱うデータフォーマットは一定時間間隔でセンサの値を取得した単変量の時系列データであることが多いです。
その際、時系列データに対してRMS(Root Mean Square)や時間微分、フーリエ変換などの様々な処理を施すことでデータの特徴や構造を把握し、効果的な前処理やアルゴリズムの検討に役立てるということがあるかと思います。
そこで、本記事ではデータの特徴を把握する上で有効な手法の一つである動的モード分解(DMD:Dynamic Mode Decomposition)の単変量時系列データへの適用方法と結果の見方についてご紹介します。

DMDは時々刻々と変化していく状態量の中に共通して存在するパターンを抽出し、そのパターンが時間にしたがってどのように振動、減衰あるいは増幅していくのかといった特徴を把握することが可能です。例えば、乱流の解析や人間の動作の解析といったところで使われています。

私自身DMDは流体の分野でよく用いられるように多次元データに適用するイメージが強かったのですが、実はセンサデータのような一次元の単変量時系列データへも適用可能だということを知ったので、今回はその点にフォーカスしてPythonでの実装と合わせてご紹介していきたいと思います。

# 実装
必要なライブラリをimportします。

```python
import numpy as np
from numpy import pi, exp, log
from numpy import dot, multiply, diag, power
from numpy.linalg import inv, eig, pinv
import matplotlib.pyplot as plt
```

## 時系列データ作成
以下の時系列データ$F$を作成します。

```python
output_dir = "analysis/"
N = 200
OBSERVASION_TIME = 1
t_raw = np.linspace(0, OBSERVASION_TIME, N)
dt = t_raw[2] - t_raw[1]
trend = 10 * (t_raw - OBSERVASION_TIME) ** 2
periodic1 = np.sin(10 * 2 * pi * t_raw) / exp(-2 * t_raw)
periodic2 = np.sin(20 * 2 * pi * t_raw)
np.random.seed(123)
noise = 1.5 * (np.random.rand(N) - 0.5)
F = trend + periodic1 + periodic2 + noise
```

$F$は以下の4つのデータの和で構築されています。
1. 緩やかに減少する2次関数
2. 10Hzで振動し減衰するsin波
3. 20Hzで振動し減衰も増幅もしないsin波
4. 一様乱数で生成されるノイズ

ここで、1をtrend、2をperiodic#1、3をperiodic#2と呼び、DMDで分解できるか検証していきます。
$F$と1~4のデータをプロットすると以下のようになります。

<figure align="center">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/befc3ed2-179b-2821-4b4f-a59b277abf4c.png" 
 alt="toy_timeseries" width="800">
</figure>

## 前処理(Hanakel行列化)
時系列データ$F=\lbrace f_0,f_1,f_2,...,f_{N-1}\rbrace$をスナップショット$\lbrace f_i, f_{i+1},...,f_{i+L-1}\rbrace\; (i=0\,,...,\,N-L)\,$に分割し、各スナップショットを連結することで以下の行列を得ます。

<figure align="center">
  <img src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/55fcb3b0-d177-4098-fb7d-c04004dd860e.png  style="text-align: center" width="600">
</figure>

前処理として時系列データを行列$D$に変換します。
これはHankel Matrixと呼ばれ、この形にデータを変形することにより単変量時系列データへのDMD適用が可能になります。Dを計算する部分の実装は参考文献[1]から引用しています。

```python
L = 70
K = N - L + 1
D = np.column_stack([F[i : i + L] for i in range(0, K)])
```

$N$はデータ数、$L$はwindow lengthと呼ばれ$2\leq L \leq N/2$で設定します。
$L$は列の長さに対応し、$K=N-L+1$が行の長さに対応します。

Dを可視化すると以下のようになります。

<figure align="center">
  <img src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/1f349b2e-e44e-edf8-8135-bfa9034c457b.png style="center" alt="Fig.2 Hankel行列" width="500">
</figure>

## DMD実行
DMDでは以下のように状態$X$から1ステップ先に進んだ状態$Y$をある線形変換$A$で表現します。
$$Y = AX \tag{1}$$
$A$の固有値解析を行うことでデータの中に存在する様々なモードの情報がわかってきます。

状態$X$と$Y$を作ります。

```python
X = D[:, :-1]
Y = D[:, 1:]
t = np.arange(0, K - 1)
```

次に$A$を求めていきます。
$A$は以下のノルムを最小化することで求めることが可能です。
$$ min\|AX-Y\|_F^2 \tag{2}$$
$\| \cdot \|_F$はフロベニウスノルムです。
上記の最小化問題の解は以下になることが知られています。
$$ A=YX^\dagger \tag{3}$$
$X^\dagger$は擬似逆行列(Moore-Penrose pseudoinverse)と呼ばれるものです。
ここで、$X$の特異値分解(SVD)を実行します。
SVDにより分解されたモードの内、重要度の低いモードを切り捨てることで効率的に情報を残したまま次元削減が可能です。

$$ X=U\Sigma V \tag{4}$$

```python
U2, Sig2, Vh2 = np.linalg.svd(X)
```

ここで各成分の寄与率累積寄与率を計算します。
$$ Relative\:contributions = \frac{\sigma_i^2}{\sum_{k=0}^{d-1}\sigma_k^2} \tag{5}$$

```python
sigma_sumsq = (Sig2 ** 2).sum()
cumlative_contrib = (Sig2 ** 2).cumsum() / sigma_sumsq * 100
```

寄与率と累積寄与率をプロットします。

<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/0e500c54-7e7c-f473-e631-4d0168d2dbdb.png" alt="eigen_plot" style="center" width="1000">
</figure>

上位5つの成分で全体の99%以上の寄与率に達しています。
今回はこれら5つの成分を使用し残りは切り捨てることにします。
ここで累積寄与率何%までの成分使うのかは特にこれといった指標はないのでデータの特徴や最終的に達成したい目的など包括的に考えて選ぶ必要があります。

以上より$A$の低次元空間での表現$\tilde{A}$は以下のようになります。
$$ \tilde{A}=U_r^* A U_r \tag{6}$$
$r$は選択した次元数なので今回は5になります。$U^*$は$U$の随伴行列です。
$V, \Sigma$は直行行列であるため、
$$ \tilde{A}= U_r^*AU_r=U_r^*AU_r\Sigma_rV_r^*V_r\Sigma_r^{-1}=U_r^*AXV\Sigma_r^{-1}=U_r^*YV_r\Sigma_r^{-1} \tag{7}$$
より$\tilde{A}$が求まります。

```python
Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig))
```
次に$\tilde{A}$の固有値と固有ベクトルを求めます。以下は$\tilde{A}$の固有値分解の式です。
$$ \tilde{A}W=W\Lambda \tag{8}$$

```python
mu, W = eig(Atil)
```

$\tilde{A}$の固有値をプロットしてみます。

<figure align="center">
  <img src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/1abc6a71-2f6f-9091-5dd4-5e9f80ad9432.png style="center" width="1000">
</figure>

各固有値の偏角はモードの振動数、絶対値はモードの減衰/増幅率に対応しています。固有値が虚数軸に対して対称になっているのは、共役な複素固有値が存在するためであり、共役な複素固有値は同じ振動特性を示すため、ここには実質3つの振動モードが存在していることがわかります。
また、各固有値の色は各固有値に対応するモードの振幅に対応させています。振幅によりそのモードがシステムの特性にどれだけ影響を与えているかを評価することが可能です。説明の都合上振幅の求め方は後述します。<br>
各固有値$\lambda_j({j=0, 1, ..., r})$の絶対値は以下のようにモードの減衰及び増幅を意味します。
* $|\lambda_j|<1:$このモードは減衰します。
* $|\lambda_j|=1:$このモードは減衰も増幅もしません。
* $|\lambda_j|>1:$このモードは増幅します。

固有値が存在している領域を拡大した右側の図を見ると以下のことが言えます。
* $\lambda_0(mode0):$実軸上に乗っているため非振動モードであり、$|\lambda_j|<1$であるため減衰します。toy dataのtrendに対応すると考えられます。
* $\lambda_1(mode1):$振動モードであり$|\lambda_j|>1$であるため増幅します。periodic#1に対応すると考えられます。
* $\lambda_3(mode3):$振動モードであり$|\lambda_j|=1$であるため減衰も増幅もしません。periodic#2に対応すると考えられます。

次に各固有値に対応する振動モードの周波数を求めます。
周波数は固有値の偏角$(\arg)$を用いて、以下のように計算されます。$\Delta t$は時系列データのタイムステップです。
$$f_j=\frac{\arg(\lambda_j)}{2\pi \Delta t} \tag{9}$$

以下のコードで各モードの振動数を出力します。

```python
freqs = np.angle(mu) / (2 * pi * dt)
for idx, freq in enumerate(freqs):
    print("mode" + str(idx) + " freq=" + str(freq))
```
出力結果

```
mode0 freq=0.0
mode1 freq=9.989363517060186
mode2 freq=-9.989363517060186
mode3 freq=19.957930838631462
mode4 freq=-19.957930838631462
```
mode0は振動数がゼロであり、trendに対応することがわかります。
periodic#1は10Hz、periodic#2は20Hzのsin波を用いて作成しており、mode1とmode2がperiodic#1、mode3とmode4がperiodic#2に一致していることがわかります。

次にDMDモードとその時間発展を計算していきます。
各モードの時間発展は以下で計算されます。
$$x(t)\approx\sum_{j=0}^{r}\phi_j\mathrm{exp}(\omega_jt)b_j=\mathbf{\Phi}\mathrm{exp}(\Omega t)\mathbf{b} \tag{10}$$

$\omega_j, \phi_j$は$A$の固有値、固有ベクトルを意味しておりそれぞれ以下で計算されます。詳しい計算過程は参考文献[2]を参照ください。
$$\mathbf{\Omega}=\mathrm{diag}(\omega_j), \quad\omega_j=\frac{\ln{(\lambda_j)}}{\Delta t} \tag{11}$$
$$ \mathbf{\Phi}=\mathbf{Y}\mathbf{V}_r\mathbf{\Sigma}_r^{-1}\mathbf{W} \tag{12}$$
$b_j$は$t=0$時の$\phi_j$方向の初期値を示して、$x(0)=X$より$X=\Phi\mathbf{b}$であるため、
$$\mathbf{b}=\mathbf{\Phi}^\dagger\mathbf{X} \tag{13}$$
で求められます。

$b_j$はDMDにおける振幅として扱われます。$b$の絶対値を先ほどの固有値プロットの色に対応させています。

ここで、(10)式の$\phi_j$がDMDモードであり、全てのスナップショット(列方向)に一貫して存在する構造になります。
(10)式の$\mathrm{exp}(\omega_j t)b_j$の部分はDMDの各モードに対応する時間発展$\Psi$でありtime dynamicsとも呼ばれます。

```python
Phi = dot(dot(dot(Y, V), inv(Sig)), W)
Psi = np.zeros([r, len(k)], dtype="complex")
for idx, mu_elem in enumerate(mu):
    for _k in k:
        Psi[idx, _k] = exp(log(mu_elem) / dt * _k * dt) * b[idx]
```

各固有値に対応するDMDモードとtime dynamicsをプロットすると以下のようになります。

<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/d5b7956b-fe23-51c0-12cc-36aea1e86bfe.png" alt="eigen_plot" style="center" width="700">
</figure>
<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/9721b55d-cdde-0ee9-cf85-919675939bce.png" alt="eigen_plot" style="center" width="700">
</figure>
<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/dc04dcca-edc7-a7f9-9260-d766ce3267a6.png" alt="eigen_plot" style="center" width="700">
</figure>
<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/0dc041db-b123-d8da-8e4d-c41b67ecd975.png" alt="eigen_plot" style="center" width="700">
</figure>
<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/b8a7226b-6d03-a476-0cf6-cc3165b05ce6.png" alt="eigen_plot" style="center" width="700">
</figure>

DMDモードはHankel行列の列方向に存在するパターンであり、DMDモードがどのように時間発展していくかを示したものがtime dynamicsとなります。以下のようなイメージを持っていただければよいかと思います。

<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/9d18ac93-1c59-9086-8e1d-92b37e58c701.png" alt="mode_dynamics_explanation" style="center" width="600">
</figure>

モード0は減衰、モード1と3は増幅、モード3と4は減衰及び増幅なしの振動であることがtime dynamicsから見て取れます。
今回のケースではHankel行列の列方向も行方向と同様に単変量(スカラー)の時間発展を意味するため、DMDモードはtime dynamics類似した傾向を示しています。
単変量以外のデータを扱う場合、ある時刻におけるシステムの状態量を1次元化したものをスナップショットとして結合したものを状態行列として扱うため、モードとDynamicsの意味するところは異なり両者を見比べながらシステムを把握することが重要になります。
ちなみに、図からも明らかなようにmode1とmode2、mode3とmode4のような共役な複素固有値は全く同じ振動特性を示すため、和を取って1つのモードと見なして問題ありません。

## 時系列データ再構築
最後に(10)式より元の時系列データの再構築を行います。

```python
# Compute DMD reconstruction for each mode
x_t_list = []
for idx, mu_elem in enumerate(mu):
    x_t = []
    for t_ in t:
        x_t.append(Phi[:, idx] * exp(log(mu_elem) / dt * t_ * dt) * b[idx])
    x_t = np.array(x_t).T
    x_t_list.append(x_t)

# Convert Hankel matrix to time-series
x_t_recon = []
for x_t in x_t_list:
    x_t_recon.append(X_to_TS(x_t))
x_t_recon = np.array(x_t_recon)
F_recon = x_t_recon.sum(axis=0)
```

元データは各モードの時間発展の線形和で求めています。X_to_TSという関数はHankel行列から時系列データに戻す関数であり参考文献[1]から引用しています。各モードの時間発展を見る必要がない場合、$\Phi$と$\Psi$の積を取るだけで時系列データ全体の再構築を計算することも可能です。
元データと各モードの再構築結果を以下に示します。

<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/4e8ee42c-2981-2cd1-8026-a6ed4224ac5f.png" alt="reconstructed_timeseries" style="center" width="800">
</figure>

適切にモード分解できていることがわかります。Periodic#1とPeriodic#2はそれぞれ2つのモードに分けられているため、一見振幅が小さくなったように見えていますが、和を取れば元の振幅に一致します。
最後に元データと再構築データを重ね合わせると以下のように非常によく一致していることがわかります。

<figure align="center">
  <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1780135/5f0c4087-5691-a8de-bee8-579498400961.png" alt="comparison" style="center", width="800">
</figure>

以上より、DMDにより単変量時系列データをモード分解し、周波数や減衰/増幅などの特徴を把握することができました。

# おわりに
このように減衰や増幅を伴う波形の解析を行う際にDMDを使用することできれいにモード分解することが可能です。
周波数成分に分解するという点でフーリエ変換と類似していますが、フーリエ変換では減衰/増幅という動的な特徴をわかりやすい形で把握することは困難です。DMDでは減衰/増幅を固有値というわかりやすい形で捉えられるため、システムを直感的に把握するという点でより有効であると考えられます。フーリエ変換とDMDの比較については、いくつかの実験を行った結果を別の記事で紹介できればと考えています。
また、今回はシンプルなトイデータを用いて検証したため、きれいな結果が出ていますが、実データに適用するときれいにモード分解されないという状況も起こりうるかと思います。実際に様々な状況でDMDによるモード分解がうまくいかないことは知られており、それらの問題に対応するためにDMDの改良版もいくつか開発されています。
今後はそういったDMDの改良版の検証や機械学習との絡め方についての検討も進めていきたいと思います。面白い検証結果が得られたらこちらも記事にしていければと考えています。

# 参考文献
[1] [Introducing SSA for Time Series Decomposition](https://www.kaggle.com/jdarcy/introducing-ssa-for-time-series-decomposition)
[2] [Dynamic Mode Decomposition in Python](http://www.pyrunner.com/weblog/2016/07/25/dmd-python/)
[3] [On Dynamic Mode Decomposition: Theory and Applications](https://arxiv.org/abs/1312.0041)
[4] [A Dynamic Mode Decomposition Approach With Hankel Blocks to Forecast Multi-Channel Temporal Series](https://ieeexplore.ieee.org/document/8718360)
[5] [Dynamic Mode Decomposition of Fast Pressure Sensitive Paint Data](https://www.mdpi.com/1424-8220/16/6/862)
[6] [Dynamic Mode Decomposition: Theory and Data Reconstruction](https://arxiv.org/abs/1909.10466)
[7] [Visualization and Selection of Dynamic Mode Decomposition Components for Unsteady Flow](https://arxiv.org/abs/2012.09633)

<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
<script type="text/x-mathjax-config">
  MathJax.Hub.Config({ tex2jax: {inlineMath: [['$', '$']]}, messageStyle: "none" });
</script>