動的モード分解にL2正則化を適用する

はじめに

$[x_0, x_1, x_2, \cdots, x_T$ ] という時系列データがあったとき、 $x_{t+1} = Ax_t$ と線形に時間発展するとします。

この線形作用素 $A$ をいくつかの固有値（時間成分）と固有モード（空間成分）の組に分解する手法が動的モード分解(Dynamic Mode Decomposition,DMD)です。

固有値の値によって、そのモードが時間的に振動するのか、増幅・減衰するのかが分かります。

僕が実装しているDMDでは、固有値や固有モードには興味がなく、線形作用素 $A$ による未来の予測に興味があります。

未来の予測を安定化させるために、DMDにL2正則化を適用しました。

普通に線形作用素 $A$ を求める

固有値や固有モードを求めたい場合は特異値分解などを用いた計算が必要ですが、 $A$ だけを求めたい場合は簡単です。

時系列データ $[x_0,x_1,x_2,\cdots,x_T$ ]があったとき、

$X_0=[x_0,x_1,\cdots,x_{T-1}$ ]

$X_1=[x_1,x_2,\cdots,x_T$ ]

とします。

ここで $x_t$ を $n$ 次元ベクトルとすると、 $X_0,X_1$ は $n \times T$ の行列になります。

線形作用素 $A$ は $X_0,X_1$ を使って

$X_1 = A X_0 + \epsilon$

という式で表すことができます。

この誤差 $\epsilon$ が最小になるように最小化問題、

$\min || X_1 - A X_0 || _2 ^2$

を解くことで $A$ を求めることができます。

この最小化問題は

$\frac{\partial L}{\partial A} = 0$

と書き表されます。（ $L = || X_1 - A X_0 || _2 ^2$ ）

$L = \mathrm{Tr} \left( \left(X_1 - A X_0 \right) ^T \left(X_1 - A X_0 \right) \right)$

$\frac{\partial L}{\partial A} = \frac{\partial}{\partial A} \mathrm{Tr} \left( \left(X_1 - A X_0 \right) ^T \left(X_1 - A X_0 \right) \right)$

$\frac{\partial L}{\partial A} = -2 \left( X_1 - A X_0 \right) X_0 ^T$

(参考:「ベクトルで微分・行列で微分」公式まとめ - Qiita )

$-2 \left( X_1 - A X_0 \right) X_0 ^T = 0$

$X_1 X_0 ^T - A X_0 X_0 ^T = 0$

$X_1 X_0 ^T = A X_0 X_0 ^T$

と式変換され、結果

$A = X_1 X_0 ^T \left( X_0 X_0 ^T \right)^{-1}$

と $A$ が求まります。

$A$ を使った予測

$A$ を使うことで未来の状態を予測できます。

時刻 $t$ の状態 $x_t$ を求めたいとき、初期状態を $x_0$ とすると

$x_t = A^{t} x_0$

と求まります。

このとき、 $A$ の固有値全てのノルムが1未満であるとき、 $x_t$ は安定となります。

(参考: 離散的な線形システムの安定性を確認する方法：制御工学入門｜Tajima Robotics)

また、理論的な事は理解してませんが、 $A$ を求める際にL2正則化を加えることで固有値の絶対値が小さくなり、安定な $A$ を求めることができるようになります。

L2正則化項を加えて $A$ を求める

$A$ を求めるための最小化問題にL2正則化項を加えると、

$\min {{|| X _ 1 - A X _ 0 || _ {2} ^ {2} } + { \lambda || A|| _ {2} ^ {2} }}$

となります。( $\lambda$ は正則化の強さを決めるパラメータ)

正則化項がない場合と同様に、

$\frac{\partial}{\partial A} \left( {|| X _ 1 - A X _ 0 || _ {2} ^ {2} } + { \lambda || A|| _ {2} ^ {2} } \right) = 0$

と、 $A$ による偏微分が0になるように求められます。

上式左辺のカッコ内第一項目の偏微分はさきほど求めたので、第二項目の偏微分

$\frac{\partial}{\partial A} { \lambda || A|| _ {2} ^ {2} }$

を求めます。

$\frac{\partial}{\partial A} { \lambda || A|| _ {2} ^ {2} } = \frac{\partial}{\partial A} { \lambda \mathrm{Tr} \left( A ^ T A \right) }$

$\mathrm{Tr} \left( A ^ T A \right) = \sum _ {i,j} a _ {ij} ^ 2$ より、

$\frac{\partial}{\partial a _ {ij}} {|| A|| _ {2} ^ {2} } = { 2 a _ {ij} }$

よって、

$\frac{\partial}{\partial A} { \lambda || A|| _ {2} ^ {2} } = { 2 \lambda A }$

となり、第一項目の偏微分と合わせると、

$\frac{\partial L}{\partial A} = -2 \left( X_1 - A X_0 \right) X_0 ^T + { 2 \lambda A } = 0$

となる。

これを解くと、

$-2 \left( X_1 - A X_0 \right) X_0 ^T + { 2 \lambda A } = 0$

$-X_1 X_0 ^T + A X_0 X_0 ^T + \lambda A = 0$

$A \left( X_0 X_0 ^T + \lambda I \right) = X_1 X_0 ^T$

$A = X_1 X_0 ^T \left( X_0 X_0 ^T + \lambda I \right) ^ {-1}$

と $A$ が求まりました。

これがL2正則化を適用した線形作用素 $A$ です。

$\lambda$ を0にすると正則化がない場合と一致することも確認できます。

結果

手持ちのデータ(100次元)でやってみました。

制御項も加わり、 $x _ {t+1} = A x _ t + B u _ t$ というような式で予測してますが結果に大きな変化はないと思うので気にしないでください。

$\lambda = 0$ の場合

$A$ のノルム(フロベニウスノルム)は394.96

この $A$ で予測をすると発散してしまいました。

f:id:marc326marcy:20200403000756p:plain — Aのヒートマップ

f:id:marc326marcy:20200403001040p:plain — 固有値の分布(横軸:実軸、縦軸:虚軸)

$\lambda = 0.1$ の場合

ノルムは49.24

これも発散してしまいました。

f:id:marc326marcy:20200403001749p:plain — Aのヒートマップ

f:id:marc326marcy:20200403001756p:plain — 固有値の分布(横軸:実軸、縦軸:虚軸)

$\lambda = 1$ の場合

ノルムは4.73

発散せず予測することができました。

f:id:marc326marcy:20200403002740p:plain — Aのヒートマップ

f:id:marc326marcy:20200403002818p:plain — 固有値の分布(横軸:実軸、縦軸:虚軸)

$\lambda = 10$ の場合

ノルムは4.38と、 $\lambda = 1$ と比較してあまり小さくなりませんでした。

$\lambda = 1$ と同じく発散しませんでした。

f:id:marc326marcy:20200403003812p:plain — Aのヒートマップ

f:id:marc326marcy:20200403003831p:plain — 固有値の分布(横軸:実軸、縦軸:虚軸)

コード

ノルムや上の画像を表示させるためのPythonコード

import numpy.linalg as la
import matplotlib.pyplot as plt

# ノルム
norm = la.norm(A)

# ヒートマップの描画
plt.imshow(A)

# 固有値取得
w, v = la.eig(A)

# 固有値の複素平面上の分布
plt.scatter(w.real, w.imag)

おわりに

おしまい

marcy's Lab

理系大学生のブログ

動的モード分解にL2正則化を適用する