動的モード分解におけるモード数低減手法の実験的比較

(1)

動的モード分解におけるモード数低減手法の実験的比較

A comparative study on sparse mode selection methods

in dynamic mode decomposition

草場彰

1∗

_{久保山哲二}

1

Akira Kusaba

1

Tetsuji Kuboyama

1

_{学習院大学計算機センター}

1

_{Computer Centre, Gakushuin University}

Abstract: Dynamic mode decomposition (DMD) is a data-driven approach for extracting oscil-latory spacial patterns from spatio-temporal data. The extracted patterns can be very useful to have an insight into the underlying mechanism of a phenomenon observed as data. However, the interpretation of the patterns becomes more diﬃcult as the number of patterns becomes larger. In this paper, we focus on two approaches for sparse mode selection to obtain a small number of pat-terns, which are based on singular value decomposition (SVD) or L1-regularization. We compared these approaches by applying them to a real-world data set of physical phenomenon to investigate which approaches are more suitable.

1 はじめに

時空間データである流体流れの解析手法として知られる「動的モード分解（Dynamic Mode Decomposition; DMD）」は時間的に繰り返す空間パターンのセットをデータ駆動で見出し，それらで流れ場を展開する手法である．各空間パターンは一定の振動数と増幅／減衰率をもって時間変化する．現象を反映していない空間パターン（正弦波基底）で展開する Fourier 解析と比較して，DMD には物理的解釈を行い易いという利点がある．しかしながら，時空間データから抽出される空間パターン数が多くなるにつれて，データの物理的解釈は難しくなる．空間パターン数は時空間データ行列のランク数だけあり，解析対象によっては，膨大な数になり得る．空間パターン数を低減するためには，特異値分解に基づいて低次元空間へデータを射影する手法と，L1 正則化により少数の空間パターンを選択する手法がある．本研究では，これら２つの手法を実際の物理現象データに適用し，実験的な比較を行った．DMD は 2008 年に Schmid と Sesterhenn によって提案されて以来 [1, 2]，流体流れ解析の研究領域 [3, 4] を越えて，神経科学 [5]，画像処理 [6, 7]，電力システム解析 [8]，スポーツ科学 [9] 等の様々な分野において活用されている．本研究では，プラズマ乱流実験において計測された電流値の時空間データ [10] を解析対象とした． ∗_{連絡先：学習院大学計算機センター} 〒 171-8588 東京都豊島区目白 1-5-1 E-mail: [email protected]

2 動的モード分解

解析対象の時空間データから，データ行列は次のように定義される． Ψ0:= [ ψ0 ψ1 . . . ψN₋₁ ] ∈ CM×N ₍₁₎ Ψ1:= [ ψ1 ψ2 . . . ψN ] ∈ CM_×N (2) ここで，ψt∈ CM は時刻 t における観測ベクトルであ り，M は観測点の数，N は時間ステップの数である． 時不変な線形システムを仮定すると，次式が成り立つ． ψt+1≈ Aψt (t = 0, . . . , N− 1) (3) 観測された時系列に最も合う行列 A は，Ψ1− AΨ0の Frobenius ノルムを最小にするものとして，次式で計算できる． A := Ψ1Ψ+0 (4) ここで，Ψ+ 0 は Ψ0の擬似逆行列であり，特異値分解 Ψ0 = U ΣV∗から，Ψ+0 = V Σ−1U∗ と計算される． 時刻 t における観測値は，初期観測値および行列 A を 用いた計算値で近似することができる．さらに行列 A の固有分解から，次のように展開される． ψt≈ Atψ0= M ∑ i=1 ϕiµtiαi (5) 人工知能学会研究会資料 SIG-FPAI-B901-01

(2)

ここで，ϕiと µiは行列 A の固有ベクトルと固有値である．それぞれ，DMD モードおよび DMD 固有値と呼ばれ，空間パターンおよび時間発展を表している．特に複素固有値の場合，対応するモードは振動特性を示 す．αiは初期観測値を反映した重ね合わせの係数（寄与率）である．データ再現の起点として初期観測値を使うことに意味はないので，最終的に寄与率は，観測された時系列と計算された時系列との差の Frobenius ノルムを最小化するように決定される．

2.1 モード数低減手法 1

ランクが r のデータ行列 Ψ0の特異値分解から求ま る行列 U ∈ CM×r_{による次元削減近似を考える．} F := U∗AU = U∗Ψ1V Σ−1 (6) の固有分解から固有ベクトル yiと固有値 µiを求めて， DMD モードを ϕi ≈ Uyi (7) と近似する．観測値は r 個の DMD モードの重ね合わ せで， ψt≈ r ∑ i=1 ϕiµtiαi (8) と表現される．十分小さな特異値をゼロとみなして，低ランク近似を行うことで，モード数の低減を図ることができる．以降，この手法を Truncated-DMD と呼ぶ．

2.2 モード数低減手法 2

Sparsity-Promoting Dynamic Mode Decomposition (SP-DMD) [11] は，式 (8) において低ランク近似を行 わず，L1 正則化によって r 個の DMD モードから重要 度の高い少数のモードを選択する手法である．観測された時系列と計算された時系列との差の Frobenius ノ ルムを最小化するための目的関数を J(α) と表記する． ここで α = (α1. . . αr) と表した．解くべき問題は次の最適化問題である． minimize α J (α) + γ r ∑ i=1 |αi| (9)

文献 [11] では，Alternating Direction Method of Mul-tipliers (ADMM) による最適化手順が示されている．最適化の結果，非ゼロ寄与率をもったモードが選択されたモードである．最終的な寄与率は，選択されなかったモードの寄与率をゼロとする拘束条件の下で，再び 正則化項を外した J(α) を最小化して決定する．

3 時空間データ

熱核融合を行うには 1 億度以上の温度が必要であり，そのような高温環境下において物質はプラズマ状態をとっている．現状では，イオン温度と電子温度をともに 1 億度まで加熱することは難しく，乱流によるエネルギーの漏れが原因と考えられている．このような背景から，核融合プラズマの乱流特性が研究されている．本研究では，直線型プラズマ生成装置において計測されたイオン飽和電流の時空間データ [10] を解析した．計測点は周方向に 64 点あり，図 1 に示すように，1 番と 64 番の計測点は 1 周して隣り合っている．サンプリングは 1 µs 間隔で行われ，410 ms にわたって計測されている．DMD による解析区間は 1 ms と設定した．図 2 は本研究で解析した或る 1 ms 区間のデータを，ヒートマップ可視化したものである．したがって，データ行 列のサイズは M = 64，N = 999 となる． 図 1: 直線型プラズマと計測システムの概略図

0.0

0.5

1.0 Time (ms)

20

40

60 Observation point

0.05

0.10

0.15

図 2: イオン飽和電流 (A) の時空間データ

(3)

4 実験

Truncated-DMD と SP-DMD によるモード数低減を 比較する．はじめに，データ行列 Ψ0の特異値を図 3 に示す．15 番目の特異値は 1 番目の特異値の 1%未満 であるので，Truncated-DMD においては，r≈ 15 と 低ランク近似を行う．一方 SP-DMD では，r = 64 で あり，パラメータ γ によって選択されるモード数が変 化する．図 4 は γ に対するモード数および誤差の程度 を表す最適化された目的関数値を示している．γ > 4 で は誤差が急激に増加していることがわかる．したがっ て，γ = 2 と設定した．このときの選択されたモード 数は 15 である． 0 10 20 30 40 50 60 103 101 101 Singular value 図 3: データ行列の特異値図 4: SP-DMD におけるモード数と誤差のトレードオフ図 5 は両手法における DMD 固有値を示している． DMD 固有値の絶対値はモードの増幅／減衰率に対応しており，1 より大きくなるほど増幅率が増大し，1 より小さくなるほど減衰率が増大する．SP-DMD によって減衰率の小さな，単位円付近の固有値が選択されていることがわかる．一方で，DMD 固有値の偏角はモードの振動数に対応している．複素共役な固有値のペアは同じ振動数のモードであり，実軸上の固有値は非振動モードである．したがって，観測データは 1 つの非振動モードと 7 つの振動モードに分解されている．Truncated-DMD 固有値は SP-つの振動モードに分解されている．Truncated-DMD 固有値と同様の領域に分布しているが，絶対値は比較的大きい．このプロットから両手法間のモードの対応関係をとることは難しい．そこで，図 6 のようなスペクトル表現を用いる．振動数において定量的な一致は見られないが，寄与率の増減に関して両スペクトルは定性的には同様の傾向を示している．したがって，両手法間のモードの対応関係をとることが可能になった． 0.2 0.4 0.6 0.8 1.0 Re 0.1 0.0 0.1 Im SP-DMD (not selected) SP-DMD (selected) 0.97 0.98 0.99 1.00 1.01 Re 0.05 0.00 0.05 Im Truncated-DMD SP-DMD (selected) 図 5: DMD 固有値 0 2 4 6 8 10 Frequency (kHz) 103 101 101 Amplitude Truncated-DMD SP-DMD (selected) 図 6: DMD スペクトル図 7 は対応すると考えられる Truncated-DMD と SP-DMD のモードごとに時間発展を比較したものである． mode-2, 4, 5 の時空間パターンは両手法間で異なって見えるが，mode-2, 5 は減衰の強いモードであり，影響は小さいと考えられる．減衰の弱いモードを見ると，両手法で同様な時空間パターンに分解されていることがわかる．また，再構成された時空間パターンも両手法間で良い一致が見られている．観測データ行列と再現データ行列の間で二乗平均平方根誤差（RMSE）を計算したところ，Truncated-DMD では RSME = 0.020 (A)， SP-DMD では RSME = 0.018 (A) であった．つまり，モードの時空間パターンから定性的に現象の物理的解

(4)

釈を行うことを目的とした場合には，どちらの手法を採用しても変わりないが，より定量的に現象のモデリングを行うことを目的とした場合には，SP-DMD の方が適していることがわかった．図 7: モードごとの時間発展（左：Truncated-DMD，右：SP-DMD）

5 むすび

本研究では，実際の物理現象データの具体例としてプラズマ乱流データを取り上げて，Truncated-DMD と SP-DMD を適用し比較した．どちらの手法においても物理的解釈が可能な範囲の少ないモード数で観測データを再構成することができた．また，僅かながら SP-DMD の方がより正確に再構成できることがわかった．一方で Truncated-DMD には，解析時間区間をシフトする場合にモードを対応付けて推移を追い易いという利点がある．今後は，Truncated-DMD における解析時間区間シフトによって，より長時間のデータ解析を進めていく．

謝辞

本研究で用いたプラズマ乱流データは，九州大学応用力学研究所稲垣滋教授から提供を受けたものです．本研究は JSPS 科研費 JP19J00871, JP19K12125 の助成を受けたものです．

参考文献

[1] Schmid, P., Sesterhenn, J.: Dynamic mode decomposition of numerical and experimental data, Bulletin of the American Physical Society, Vol. 53 (2008)

[2] Schmid, P. J.: Dynamic mode decomposition of numerical and experimental data, Journal of

Fluid Mechanics, Vol. 656, pp. 5–28 (2010)

[3] Mezić, I.: Analysis of fluid flows via spectral properties of the koopman operator, Annual

Re-view of Fluid Mechanics, Vol. 45, pp. 357–378

(2013)

[4] Tu, J. H., Rowley, C. W., Luchtenburg, D. M., Brunton, S. L., Kutz, J. N.: On dynamic mode decomposition: theory and applications, arXiv

preprint arXiv: 1312.0041 (2013)

[5] Brunton, B. W., Johnson, L. A., Ojemann, J. G., Kutz, J. N.: Extracting spatial–temporal co-herent patterns in large-scale neural recordings using dynamic mode decomposition, Journal of

Neuroscience Methods, Vol. 258, pp. 1–15 (2016)

[6] Kutz, J. N., Fu, X., Brunton, S. L.: Mul-tiresolution dynamic mode decomposition, SIAM

Journal on Applied Dynamical Systems, Vol. 15,

(5)

[7] Takeishi, N., Kawahara, Y., Yairi, T.: Sparse nonnegative dynamic mode decomposition, In

2017 IEEE International Conference on Image Processing, pp. 2682–2686 (2017)

[8] Susuki, Y., Mezic, I., Raak, F., Hikihara, T.: Ap-plied koopman operator theory for power systems technology, Nonlinear Theory and Its

Applica-tions, IEICE, Vol. 7, pp. 430–459 (2016)

[9] Fujii, K., Inaba, Y., Kawahara, Y.: Koopman spectral kernels for comparing complex dynam-ics: Application to multiagent sport plays, In

Joint European Conference on Machine Learn-ing and Knowledge Discovery in Databases, pp.

127–139 (2017)

[10] Yamada, T., Itoh, S. I., Maruta, T., Kasuya, N., Nagashima, Y., Shinohara, S., Terasaka, K., Yagi, M., Inagaki, S., Kawai, Y., et al.: Anatomy of plasma turbulence, Nature Physics, Vol. 4, pp. 721–725 (2008)

[11] Jovanovi´c, M. R., Schmid, P. J., Nichols, J. W.: Sparsity-promoting dynamic mode decom-position, Physics of Fluids, Vol. 26, p. 024103 (2014)

動的モード分解におけるモード数低減手法の実験的比較