はじめに

前回投稿のDynamic Mode Decompositionの続き。
本稿ではSINDyの導入に関する概要を理解し、ローレンツ系を用いてシステムの推定のデモを実施する。早く次章のクープマン作用素に行きたいので、応用や拡張、歴史的背景はすっ飛ばす。

概要

SINDy（Sparse Identification of Nonlinear Dynamics）とは、物理システムのダイナミクスを支配する非線形な微分方程式を、データから自動的に抽出するアルゴリズムである。ただし、ハイクオリティな測定値とデータ数が求められる。
後のデモコードでは辞書行列の作成をスクラッチで行っているが、オープンソースのPySINDyというのがあり、メソッドが用意されているので、比較的容易に実装できる。

SINDy: Sparse Identification of Nonlinar Dynamics

非線形ダイナミクスの疎同定アルゴリズム: SINDy は、多くの動的システムが以下の特性を持つことを利用して、候補となりうる全モデル構造とその最適パラメータの組み合わせ探索を回避する。

\begin{align}
\frac{d}{d\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})
\end{align}

ここで、右辺関数$\boldsymbol{f}(\boldsymbol{x})$が有効項を少数しか持たないことに着目する。

例えば、ローレンツ系を考えたとき、

\begin{align}
\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}
\end{align}

各方程式で大きな影響を与える項（e.x.$-xz$、$xy$など）はせいぜい1~2個である。
運動の遅い系では、摩擦力（上の例だと$\sigma$）、運動が速い系だと慣性（二階微分）が支配的となる。

次に、一般化線形モデルで関数$\boldsymbol{f}$を近似する。
適当な$p$個の非線形関数$\lbrace\theta_k(\boldsymbol{x})\rbrace_{k=1}^p$と、非ゼロ係数ベクトル$\boldsymbol{\xi}=\left(\xi_1\cdots\xi_p\right)^*$を用いて、

\begin{align}
\boldsymbol{f}(\boldsymbol{x})\simeq\sum_{k=1}^p\theta_k(\boldsymbol{x})\xi_k
=
\begin{bmatrix}
| & & | \\
\theta_1(\boldsymbol{x}) & \cdots & \theta_p(\boldsymbol{x}) \\
| & & |
\end{bmatrix}
\begin{bmatrix}
\xi_1 \\
\vdots \\
\xi_p
\end{bmatrix}
=\Theta(\boldsymbol{x})\boldsymbol{\xi}
\end{align}

この定式化は力学系の項数にペナルティを与えるスパース回帰を用いて、系の中で重要な項を解くことができる。

さて、式$(1)$の微分方程式$\frac{d}{d\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$の解として得られる時系列データ$\boldsymbol{x}(t)$について考える。この時系列データを行列で表現すると、

\begin{align}
\boldsymbol{X}=
\begin{bmatrix}
\boldsymbol{x}(t_1) & \boldsymbol{x}(t_2) & \cdots & \boldsymbol{x}(t_m)
\end{bmatrix}^T
\end{align}

同様に、状態ベクトル$\boldsymbol{x}(t)$の一階微分を行列として表現

\begin{align}
\dot{\boldsymbol{X}}=
\begin{bmatrix}
\dot{\boldsymbol{x}}(t_1) & \dot{\boldsymbol{x}}(t_2) & \cdots & \dot{\boldsymbol{x}}(t_m)
\end{bmatrix}^T
\end{align}

次に、データ行列$\boldsymbol{X}$から非線形関数のセット$\boldsymbol{\Theta}(\boldsymbol{X})=\lbrace\theta_k(\boldsymbol{X})\rbrace_{k=1}^\infty$を構築する。

\begin{align}
\boldsymbol{\Theta}(\boldsymbol{X})=
\begin{bmatrix}
\boldsymbol{1} & \boldsymbol{X} & \boldsymbol{X}^2 & \cdots & \boldsymbol{X}^d & \cdots & \sin(\boldsymbol{X}) & \cdots
\end{bmatrix}
\end{align}

ここで、$\boldsymbol{\Theta}(\boldsymbol{X})$は$\boldsymbol{X}$に基づく候補となりうる全ての非線形関数を列ベクトルに持つ行列である。具体例として以下のような成分を持つ。

$\boldsymbol{1}$: 定数項
$\boldsymbol{X}$: 元のデータ行列
$\boldsymbol{X}^d$: $\boldsymbol{X}$の多項式関数
$\sin(\boldsymbol{X})$: $\boldsymbol{X}$の三角関数

すると、元の微分方程式$\frac{d}{d\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$はデータ行列を用いて以下のように表現できる。

\begin{align}
\dot{\boldsymbol{X}}&=
\boldsymbol{\Theta}(\boldsymbol{X})\boldsymbol{\Xi} \\
&=
\begin{bmatrix}
\boldsymbol{1} & \cdots & \boldsymbol{X}^d & \cdots & \sin(\boldsymbol{X}) & \cdots
\end{bmatrix}
\begin{bmatrix}
| & & | \\
\boldsymbol{\xi}_1 & \cdots & \boldsymbol{\xi}_k & \cdots \\
| & & |
\end{bmatrix}
\end{align}

この$\boldsymbol{\xi}_k$は時系列データ$\boldsymbol{x}(t)$の$k$番目のアクティブ項を決める係数ベクトルである。
倹約的なモデルでは、$\boldsymbol{\Xi}$をなるべく少ない列で$\boldsymbol{x}(t)$にフィッティングするモデルを考える。このようなモデルは、凸$\ell_1$正規化された疎回帰で識別できる。

\begin{align}
\tilde{\boldsymbol{\xi}_k}=
\mathop{\rm arg~min}\limits_{\boldsymbol{\xi}_k}\|\dot{\boldsymbol{X}}_k-\boldsymbol{\Theta}(\boldsymbol{X})\boldsymbol{\xi}_k\|_2+\lambda\|\boldsymbol{\xi}_k\|_1
\end{align}

ここで$\dot{\boldsymbol{X}}_k$は$\dot{\boldsymbol{X}}$の$k$列目、$\lambda$は正則化係数である。
この優決定性問題（$m>>n$）の数値安定性向上のためにかつては圧縮センシングが用いられていたが、ここではLASSOやSTLS（Sequential Threshold Least-Squares）を用いる。

ローレンツ系の動的モデルの推定

例として、ローレンツ系の微分方程式をSINDy（STLS）を用いて特定する。
先ずはライブラリのインポートから。

from itertools import combinations_with_replacement

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
plt.rcParams['figure.figsize'] = [12,8]
plt.rcParams.update({'font.size': 12})

続いて、ローレンツ系を適当に定義する。

# Simulate the Lorenz System
sigma = 10.0
beta = 2.0
rho = 25.0


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.0, 1.0, 20.0)
dt = 0.01
t = np.arange(0, 50, dt)

次にscipy.integrate.odeintを用いて定義したローレンツ系の微分方程式を解き、位置$\boldsymbol{X}\in\mathbb{R}^{t\times 3}$を得る。

X = odeint(lorenz, X0, t, args=(sigma, beta, rho))

さらに、この値を使って元のローレンツ系の関数に代入し、時間微分係数を計算

# Compute derivative
dX = np.zeros_like(X)
for j in range(len(t)):
    dX[j,:] = lorenz(X[j,:], 0, sigma, beta, rho)

非線形関数の辞書$\boldsymbol{\Theta}$を生成する関数を定義する。これは、探索範囲の次数${\rm polyorder}$を与えれば、$\boldsymbol{\Theta} \in \mathbb{R}^{t \times \sum\scriptscriptstyle \binom{i+2}{i}}$を返す関数である。
例えば、以下の関数でpoolData(X, 3, 3)を与えると、行方向は

\begin{matrix}
1 \\
x \\
y \\
z \\
x^2 \\
xy \\
xz \\
y^2 \\
yz \\
z^2 \\
x^3 \\
\vdots \\
z^3
\end{matrix}

のようになる。

# Generate the library of candidate nonlinear functions
def poolData(yin, nVars, polyorder):
    n = yin.shape[0]
    yout = np.ones((n, 1))  # poly order 0
    labels = ['1']  # List of labels (first constant term)

    # poly order 1 to N
    for order in range(1, polyorder + 1):
        for comb in combinations_with_replacement(range(nVars), order):
            yout = np.append(yout, np.prod(yin[:, comb], axis=1).reshape(n, 1), axis=1)
            # Create and add labels (mapping variable names to x, y, z)
            labels.append(''.join([f'x{idx+1}' for idx in comb]))

    return yout, labels

SINDyを実装。STLSを用いてスパースなモデルを見つける。
STLSについては計算手順について詳細を記す。

Step.1

初めに最小二乗法を用いて初期回帰係数$\boldsymbol{\Xi}^{(0)}$を計算する。

\begin{align}
\boldsymbol{\Xi}^{(0)}=
\mathop{\rm arg~min}\limits_{\boldsymbol{\Xi}}
\|\dot{\boldsymbol{X}}-\boldsymbol{\Theta\Xi}\|
\end{align}

Step. 2

続いて、回帰係数の値が小さいものをゼロに置換することで閾値に基づくスパース化を行う。これが$\ell_1$正則化項の役目を担う。

\begin{align}
\xi_{i,j}^{(1)}=
\begin{cases}
0,\quad&\text{if }|\xi_{i,j}^{(0)}|<\lambda \\
\xi_{i,j}^{(0)},\quad&\text{otherwise}\\
\end{cases}
\end{align}

Step.3

Step.2のスパース化の操作を施した$\boldsymbol{\Xi}^{(1)}$に対して、非ゼロの係数のみに着目して、列ごとに最小二乗法を施す。

\begin{align}
\boldsymbol{\xi}_{\boldsymbol{S},j}^{(2)}=
\mathop{\rm arg~min}\limits_{\boldsymbol{\xi}_{\boldsymbol{S},j}^{(1)}}
\|\dot{\boldsymbol{X}_{,j}}-\boldsymbol{\Theta}_{,\boldsymbol{S}}\boldsymbol{\xi}_{\boldsymbol{S},j}^{(1)}\|
\end{align}

ただし、ここで$\boldsymbol{S}$は列$j$における、非ゼロの係数インデックスを表す。

Step.4

上記、Step.2~Step.3の操作を一定の回数繰り返す。あるいは係数の収束判定を設けて、その条件を満足するまで繰り返す。

# SINDy function definitions
def sparseDynamics(Theta, dXdt, lamb, dim, trials=10):
    Xi = np.linalg.lstsq(Theta, dXdt, rcond=None)[0]  # Initial estimation

    for _ in range(trials):
        smallinds = np.abs(Xi) < lamb
        Xi[smallinds] = 0  # find small coef.
        for ind in range(dim):
            biginds = smallinds[:, ind] == False
            Xi[biginds, ind] = np.linalg.lstsq(
                Theta[:, biginds], dXdt[:, ind], rcond=None
            )[0]

    return Xi

Theta, labels = poolData(X, nVars=3, polyorder=3)
lamb = 0.025
Xi = sparseDynamics(Theta, dX, lamb, n=3)
Xi

実行結果は以下の通り。

array([[  0.,   0.,   0.],
       [-10.,  25.,   0.],
       [ 10.,  -1.,   0.],
       [  0.,   0.,  -2.],
       [  0.,   0.,   0.],
       [  0.,   0.,   1.],
       [  0.,  -1.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.],
       [  0.,   0.,   0.]])

行方向は先述の通りだが、列方向はそれぞれ$\left(\frac{dx}{dt}\ \frac{dy}{dt}\ \frac{dz}{dt}\right)$を表している。
つまり、

\begin{align}
\cfrac{dx}{dt}&=-10x+10y \\
\cfrac{dy}{dt}&=25x-y-xz \\
\cfrac{dz}{dt}&=-2z+xy
\end{align}

と、元の微分方程式を同定できたことがわかる。

おわりに

今回は急ぎのため、一旦ここで終わることにする。
後半部分を流し見したところ、偏微分方程式の探索法があるらしくこれは別のタイミングで是非読みたい。