物理現象を記述する支配方程式（微分方程式）は、これまで専門家の洞察と理論的考察によって導かれてきました。しかし、近年の計測技術の発展により大量の時系列データが取得可能になった今、「データから支配方程式そのものを発見できないか？」という問いが注目されています。
SINDy（Sparse Identification of Nonlinear Dynamics）は、まさにこの問いに答える手法です。観測データから非線形力学系の支配方程式をスパース回帰によって同定します。物理インフォームドMLの中でも、解釈可能性が高く、発見された方程式をそのまま理論的に分析できる点が大きな強みです。
本記事の内容
- SINDyの動機と基本的な考え方
- 候補関数ライブラリの構成（多項式・三角関数・指数関数）
- 数値微分の手法（有限差分・多項式フィッティング）
- スパース回帰（STLS）アルゴリズムの導出
- LASSOとの比較とノイズに対するロバスト性
- PDE-FIND への拡張の概要
- PythonによるLorenz方程式からのSINDy再構成
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
SINDyとは
SINDy（Sparse Identification of Nonlinear Dynamics）は、Bruntonら（2016）によって提案された手法で、時系列データから力学系の支配方程式を発見するフレームワークです。
多くの物理系は、状態変数 $\bm{x}(t) \in \mathbb{R}^n$ に対して次の形の常微分方程式で記述されます。
$$ \frac{d\bm{x}}{dt} = \bm{f}(\bm{x}) $$
SINDyの核心的なアイデアは、「自然界の支配方程式は、可能な関数の空間の中でスパース（疎）な表現を持つ」という仮定です。つまり、$\bm{f}(\bm{x})$ は多数の候補関数のうち、ごく少数の関数の線形結合で表現できるはずだ、ということです。
例えば、単振り子の運動方程式 $\ddot{\theta} = -\frac{g}{l}\sin\theta$ は、$\{1, \theta, \theta^2, \sin\theta, \cos\theta, \dots\}$ という候補関数の中から $\sin\theta$ の項だけが選ばれた「スパースな」表現です。
SINDyの数学的定式化
データ行列の構成
$m$ 個の時刻 $t_1, t_2, \dots, t_m$ における状態変数の観測データを行列にまとめます。
$$ \bm{X} = \begin{pmatrix} \bm{x}(t_1)^T \\ \bm{x}(t_2)^T \\ \vdots \\ \bm{x}(t_m)^T \end{pmatrix} \in \mathbb{R}^{m \times n} $$
同様に、各時刻における時間微分を行列にまとめます。
$$ \dot{\bm{X}} = \begin{pmatrix} \dot{\bm{x}}(t_1)^T \\ \dot{\bm{x}}(t_2)^T \\ \vdots \\ \dot{\bm{x}}(t_m)^T \end{pmatrix} \in \mathbb{R}^{m \times n} $$
候補関数ライブラリ行列
次に、候補関数のライブラリ行列 $\bm{\Theta}(\bm{X})$ を構成します。各列が一つの候補関数に対応します。
$$ \bm{\Theta}(\bm{X}) = \begin{pmatrix} | & | & | & | & | & \cdots \\ \bm{1} & \bm{X} & \bm{X}^{P_2} & \bm{X}^{P_3} & \sin(\bm{X}) & \cdots \\ | & | & | & | & | & \cdots \end{pmatrix} \in \mathbb{R}^{m \times p} $$
ここで $\bm{X}^{P_2}$ は2次の多項式項（$x_1^2, x_1 x_2, x_2^2, \dots$）、$\bm{X}^{P_3}$ は3次の多項式項を表します。$p$ は候補関数の総数です。
具体的に、状態変数が $\bm{x} = (x_1, x_2, x_3)^T$ の3次元系で多項式3次までと三角関数を含める場合、候補関数は以下のようになります。
$$ \bm{\Theta}(\bm{X}) = \begin{pmatrix} 1 & x_1 & x_2 & x_3 & x_1^2 & x_1 x_2 & x_1 x_3 & x_2^2 & x_2 x_3 & x_3^2 & x_1^3 & \cdots & \sin(x_1) & \cos(x_1) & \cdots \end{pmatrix} $$
スパース係数行列
SINDyの目標は、次の線形方程式を満たすスパースな係数行列 $\bm{\Xi}$ を見つけることです。
$$ \dot{\bm{X}} = \bm{\Theta}(\bm{X}) \bm{\Xi} $$
ここで $\bm{\Xi} \in \mathbb{R}^{p \times n}$ はスパース係数行列であり、$\bm{\Xi}$ の第 $k$ 列 $\bm{\xi}_k$ は、$\dot{x}_k$ がどの候補関数の組み合わせで表されるかを示します。
$$ \dot{x}_k = \bm{\Theta}(\bm{X}) \bm{\xi}_k, \quad k = 1, 2, \dots, n $$
候補関数ライブラリの設計
候補関数ライブラリの設計はSINDyの性能を大きく左右します。含めるべき候補関数の種類と、その構成方法を整理します。
多項式基底
最も基本的な候補関数は多項式です。$n$ 次元状態空間で $d$ 次までの多項式を含めると、候補関数の数は
$$ p_{\text{poly}} = \binom{n + d}{d} $$
となります。例えば $n = 3, d = 3$ の場合、定数項を含めて $\binom{6}{3} = 20$ 個の候補関数になります。
三角関数基底
振動系や周期的な現象を記述する場合、三角関数を含めることが有効です。
$$ \{\sin(x_i), \cos(x_i), \sin(2x_i), \cos(2x_i), \dots\}, \quad i = 1, \dots, n $$
指数関数基底
減衰や成長を含む系では指数関数が有用です。
$$ \{e^{x_i}, e^{-x_i}\}, \quad i = 1, \dots, n $$
候補関数の数とスパース性のトレードオフ
候補関数を増やすほど表現力は高まりますが、以下の問題が生じます。
- 計算コスト: $p$ が大きいほど回帰問題のサイズが増加
- 偽発見: 候補が多いと偶然フィットする偽の項が選ばれやすい
- 条件数の悪化: $\bm{\Theta}(\bm{X})$ の条件数が大きくなり、数値的に不安定になる
実用的には、対象系の事前知識に基づいて候補関数を選ぶことが推奨されます。
数値微分の手法
SINDyを適用するためには、観測データ $\bm{X}$ から時間微分 $\dot{\bm{X}}$ を推定する必要があります。数値微分はノイズに対して非常に敏感であり、SINDyの精度を大きく左右します。
有限差分法
最も単純な方法は有限差分です。中心差分を用いると、
$$ \dot{x}(t_i) \approx \frac{x(t_{i+1}) – x(t_{i-1})}{2\Delta t} $$
ここで $\Delta t = t_{i+1} – t_i$ は時間刻みです。中心差分の打切り誤差は $O(\Delta t^2)$ です。
4次精度の中心差分を用いるとさらに精度が向上します。
$$ \dot{x}(t_i) \approx \frac{-x(t_{i+2}) + 8x(t_{i+1}) – 8x(t_{i-1}) + x(t_{i-2})}{12\Delta t} $$
この打切り誤差は $O(\Delta t^4)$ です。
多項式フィッティング
局所的な多項式フィッティング（Savitzky-Golayフィルタ）は、ノイズの影響を軽減しながら微分を推定する手法です。
時刻 $t_i$ の周辺 $2w+1$ 点のデータに $d$ 次多項式をフィットし、その微分値を用いる方法です。
$$ p(t) = a_0 + a_1(t – t_i) + a_2(t – t_i)^2 + \cdots + a_d(t – t_i)^d $$
最小二乗法で係数 $a_0, a_1, \dots, a_d$ を求めると、$\dot{x}(t_i) \approx a_1$ が微分の推定値になります。
この方法は窓幅 $w$ をパラメータとして、平滑化の度合いを制御できます。窓幅が大きいほどノイズに対してロバストになりますが、高周波の信号成分も失われます。
全変動正則化微分（TVD）
ノイズが大きい場合に有効な方法として、全変動正則化（Total Variation Regularized Differentiation）があります。微分 $u = \dot{x}$ を次の最適化問題として定式化します。
$$ \min_u \int |u'(t)|\ dt \quad \text{subject to} \quad \int_0^t u(s)\ ds = x(t) – x(0) $$
これは $u$ の全変動を最小化しつつ、$u$ の積分が元の信号 $x(t)$ を再現するという制約です。離散化すると、
$$ \min_{\bm{u}} \|\bm{D}\bm{u}\|_1 + \frac{\mu}{2}\|\bm{A}\bm{u} – \bm{x}\|_2^2 $$
ここで $\bm{D}$ は差分行列、$\bm{A}$ は積分行列（下三角行列に $\Delta t$ を乗じたもの）、$\mu$ は正則化パラメータです。
STLS（Sequential Thresholded Least Squares）アルゴリズム
SINDyの中核となるスパース回帰アルゴリズムがSTLSです。
問題設定
各状態変数 $k$ について、次の最適化問題を解きます。
$$ \min_{\bm{\xi}_k} \|\dot{\bm{X}}_k – \bm{\Theta}(\bm{X})\bm{\xi}_k\|_2^2 \quad \text{subject to} \quad \|\bm{\xi}_k\|_0 \leq s $$
ここで $\dot{\bm{X}}_k$ は $\dot{\bm{X}}$ の第 $k$ 列、$\|\cdot\|_0$ は $\ell_0$ ノルム（非ゼロ要素の数）です。しかし $\ell_0$ 制約の最適化はNP困難であるため、STLSは反復的な閾値処理で近似解を求めます。
STLSアルゴリズムの導出
ステップ1: 通常の最小二乗法で初期解を求めます。
$$ \bm{\xi}_k^{(0)} = (\bm{\Theta}^T \bm{\Theta})^{-1} \bm{\Theta}^T \dot{\bm{X}}_k $$
これは正規方程式の解であり、導出は以下の通りです。残差の二乗和を最小化する問題
$$ J(\bm{\xi}_k) = \|\dot{\bm{X}}_k – \bm{\Theta}\bm{\xi}_k\|_2^2 = (\dot{\bm{X}}_k – \bm{\Theta}\bm{\xi}_k)^T(\dot{\bm{X}}_k – \bm{\Theta}\bm{\xi}_k) $$
を展開すると、
$$ \begin{align} J(\bm{\xi}_k) &= \dot{\bm{X}}_k^T \dot{\bm{X}}_k – 2\dot{\bm{X}}_k^T \bm{\Theta}\bm{\xi}_k + \bm{\xi}_k^T \bm{\Theta}^T \bm{\Theta} \bm{\xi}_k \end{align} $$
$\bm{\xi}_k$ で微分してゼロとおくと、
$$ \begin{align} \frac{\partial J}{\partial \bm{\xi}_k} &= -2\bm{\Theta}^T \dot{\bm{X}}_k + 2\bm{\Theta}^T \bm{\Theta} \bm{\xi}_k = \bm{0} \\ \therefore \quad \bm{\xi}_k &= (\bm{\Theta}^T \bm{\Theta})^{-1} \bm{\Theta}^T \dot{\bm{X}}_k \end{align} $$
ステップ2: 閾値 $\lambda$ より小さい係数をゼロにします。
$$ \xi_{k,j}^{(i+1)} = \begin{cases} \xi_{k,j}^{(i)} & \text{if } |\xi_{k,j}^{(i)}| \geq \lambda \\ 0 & \text{if } |\xi_{k,j}^{(i)}| < \lambda \end{cases} $$
ステップ3: 残った非ゼロ要素のインデックス集合 $\mathcal{S}^{(i+1)} = \{j : \xi_{k,j}^{(i+1)} \neq 0\}$ に対して、再度最小二乗法を適用します。
$$ \bm{\xi}_{k,\mathcal{S}^{(i+1)}}^{(i+1)} = (\bm{\Theta}_{\mathcal{S}^{(i+1)}}^T \bm{\Theta}_{\mathcal{S}^{(i+1)}})^{-1} \bm{\Theta}_{\mathcal{S}^{(i+1)}}^T \dot{\bm{X}}_k $$
ここで $\bm{\Theta}_{\mathcal{S}}$ は $\bm{\Theta}$ から $\mathcal{S}$ に対応する列のみを抽出した部分行列です。
ステップ4: ステップ2-3を収束するまで繰り返します。収束条件は、非ゼロ要素の集合が変化しなくなることです。
STLSの擬似コード
Algorithm: STLS (Sequential Thresholded Least Squares)
Input: Θ (ライブラリ行列), Ẋ_k (微分データ), λ (閾値)
Output: ξ_k (スパース係数ベクトル)
1: ξ_k ← (Θ^T Θ)^{-1} Θ^T Ẋ_k // 初期最小二乗解
2: repeat
3: small_inds ← {j : |ξ_{k,j}| < λ}
4: ξ_k[small_inds] ← 0
5: big_inds ← {j : |ξ_{k,j}| ≥ λ}
6: ξ_k[big_inds] ← (Θ_{big_inds}^T Θ_{big_inds})^{-1} Θ_{big_inds}^T Ẋ_k
7: until convergence
8: return ξ_k
LASSOとの比較
スパース回帰の代表格であるLASSO（Least Absolute Shrinkage and Selection Operator）との違いを明確にしておきましょう。
LASSOの定式化
LASSOは $\ell_1$ 正則化を用いたスパース回帰です。
$$ \min_{\bm{\xi}_k} \frac{1}{2m}\|\dot{\bm{X}}_k – \bm{\Theta}\bm{\xi}_k\|_2^2 + \alpha\|\bm{\xi}_k\|_1 $$
ここで $\alpha > 0$ は正則化パラメータ、$\|\bm{\xi}_k\|_1 = \sum_j |\xi_{k,j}|$ は $\ell_1$ ノルムです。
STLSとLASSOの比較
| 特性 | STLS | LASSO |
|---|---|---|
| 正則化 | ハード閾値処理 | $\ell_1$ ペナルティ |
| バイアス | 非ゼロ係数はバイアスなし | 係数が縮小されるバイアスあり |
| ハイパーパラメータ | 閾値 $\lambda$ | 正則化パラメータ $\alpha$ |
| 計算コスト | 反復最小二乗（軽量） | 凸最適化（やや重い） |
| 理論的保証 | 限定的 | RIPなどの理論的保証あり |
STLSの大きな利点は、非ゼロ係数にバイアスが生じないことです。LASSOでは $\ell_1$ ペナルティにより全ての係数が原点方向に縮小されるため、発見された方程式の係数が真の値からずれる傾向があります。SINDyでは物理法則の正確な係数値が重要であるため、STLSが適しています。
ノイズに対するロバスト性
実データにはノイズが含まれるため、SINDyのノイズロバスト性は実用上重要です。
ノイズの影響の数学的分析
観測データにガウスノイズ $\bm{\epsilon} \sim \mathcal{N}(\bm{0}, \sigma^2 \bm{I})$ が加わった場合、
$$ \bm{X}_{\text{noisy}} = \bm{X}_{\text{true}} + \bm{\epsilon} $$
このノイズは2つの経路でSINDyに影響します。
- 数値微分の増幅: 微分演算はノイズを $O(\sigma / \Delta t)$ に増幅します
- ライブラリ行列の汚染: $\bm{\Theta}(\bm{X}_{\text{noisy}}) \neq \bm{\Theta}(\bm{X}_{\text{true}})$
特に数値微分の影響が大きく、ノイズレベルが高い場合は前述のTVD微分やSavitzky-Golayフィルタの適用が必須になります。
ロバスト性を向上させる手法
- アンサンブルSINDy: 複数のサブサンプルに対してSINDyを適用し、安定して選ばれる項だけを採用する
- 積分形SINDy: 微分ではなく積分形式で定式化し、ノイズ増幅を回避する
$$ \bm{x}(t) – \bm{x}(0) = \int_0^t \bm{\Theta}(\bm{x}(s)) \bm{\Xi}\ ds $$
PDE-FINDへの拡張
SINDyのフレームワークは偏微分方程式（PDE）の発見にも拡張できます。これがPDE-FIND（PDE Functional Identification of Nonlinear Dynamics）です。
偏微分方程式の定式化
PDEの一般形を考えます。
$$ u_t = f(u, u_x, u_{xx}, u_{xxx}, \dots, u u_x, u^2, \dots) $$
SINDyと同様に、候補関数ライブラリには空間微分も含めます。
$$ \bm{\Theta} = \begin{pmatrix} \bm{1} & \bm{U} & \bm{U}_x & \bm{U}_{xx} & \bm{U}^2 & \bm{U}\bm{U}_x & \cdots \end{pmatrix} $$
同じSTLSアルゴリズムを適用して、
$$ \bm{U}_t = \bm{\Theta}\bm{\xi} $$
を解くことで、PDEの形を発見できます。例えば、Burgers方程式 $u_t = -uu_x + \nu u_{xx}$ やKdV方程式 $u_t = -uu_x – u_{xxx}$ がデータから正しく同定された例が報告されています。
PythonによるSINDy実装
Lorenz方程式のデータ生成とSINDy再構成
Lorenz方程式は以下の3次元非線形力学系です。
$$ \begin{align} \dot{x} &= \sigma(y – x) \\ \dot{y} &= x(\rho – z) – y \\ \dot{z} &= xy – \beta z \end{align} $$
標準的なパラメータは $\sigma = 10, \rho = 28, \beta = 8/3$ です。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from itertools import combinations_with_replacement
# ===== 1. Lorenz方程式のデータ生成 =====
def lorenz(t, state, sigma=10.0, rho=28.0, beta=8.0/3.0):
"""Lorenz方程式の右辺"""
x, y, z = state
dxdt = sigma * (y - x)
dydt = x * (rho - z) - y
dzdt = x * y - beta * z
return [dxdt, dydt, dzdt]
# 初期条件と時間設定
x0 = [-8.0, 7.0, 27.0]
t_span = (0, 10)
dt = 0.001
t_eval = np.arange(t_span[0], t_span[1], dt)
# 数値積分
sol = solve_ivp(lorenz, t_span, x0, t_eval=t_eval, method='RK45',
rtol=1e-12, atol=1e-12)
X = sol.y.T # (m, 3)
t = sol.t
print(f"データ点数: {X.shape[0]}, 状態変数の次元: {X.shape[1]}")
# Lorenzアトラクタの可視化
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.plot(X[:, 0], X[:, 1], X[:, 2], lw=0.5, alpha=0.8)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('Lorenz Attractor')
plt.tight_layout()
plt.show()
# ===== 2. 数値微分（4次精度中心差分） =====
def numerical_derivative(X, dt):
"""4次精度中心差分による数値微分"""
m, n = X.shape
dXdt = np.zeros_like(X)
# 内部点: 4次精度中心差分
for i in range(2, m - 2):
dXdt[i] = (-X[i+2] + 8*X[i+1] - 8*X[i-1] + X[i-2]) / (12 * dt)
# 境界点: 2次精度中心差分
dXdt[1] = (X[2] - X[0]) / (2 * dt)
dXdt[-2] = (X[-1] - X[-3]) / (2 * dt)
# 端点: 前進/後退差分
dXdt[0] = (-3*X[0] + 4*X[1] - X[2]) / (2 * dt)
dXdt[-1] = (3*X[-1] - 4*X[-2] + X[-3]) / (2 * dt)
return dXdt
dXdt = numerical_derivative(X, dt)
print(f"数値微分の形状: {dXdt.shape}")
# 解析的な微分値との比較（最初の点で確認）
dXdt_true = np.array([lorenz(0, X[0])])
print(f"数値微分 (t=0): {dXdt[0]}")
print(f"解析微分 (t=0): {dXdt_true[0]}")
print(f"相対誤差: {np.abs(dXdt[0] - dXdt_true[0]) / (np.abs(dXdt_true[0]) + 1e-10)}")
# ===== 3. 候補関数ライブラリの構成 =====
def build_library(X, poly_order=3, include_trig=False):
"""候補関数ライブラリ Θ(X) を構成する"""
m, n = X.shape
library = []
labels = []
# 定数項
library.append(np.ones(m))
labels.append('1')
# 1次の項
var_names = ['x', 'y', 'z']
for i in range(n):
library.append(X[:, i])
labels.append(var_names[i])
# 2次以上の多項式項
for order in range(2, poly_order + 1):
for combo in combinations_with_replacement(range(n), order):
term = np.ones(m)
label_parts = []
for idx in combo:
term *= X[:, idx]
label_parts.append(var_names[idx])
library.append(term)
labels.append(''.join(label_parts))
# 三角関数（オプション）
if include_trig:
for i in range(n):
library.append(np.sin(X[:, i]))
labels.append(f'sin({var_names[i]})')
library.append(np.cos(X[:, i]))
labels.append(f'cos({var_names[i]})')
Theta = np.column_stack(library)
return Theta, labels
Theta, labels = build_library(X, poly_order=3, include_trig=False)
print(f"ライブラリ行列の形状: {Theta.shape}")
print(f"候補関数 ({len(labels)}個): {labels}")
# ===== 4. STLS（Sequential Thresholded Least Squares）アルゴリズム =====
def stls(Theta, dXdt_col, threshold=0.1, max_iter=20):
"""
Sequential Thresholded Least Squares
Theta: ライブラリ行列 (m, p)
dXdt_col: 微分データの1列 (m,)
threshold: 閾値 λ
max_iter: 最大反復回数
"""
# ステップ1: 初期最小二乗解
xi = np.linalg.lstsq(Theta, dXdt_col, rcond=None)[0]
for iteration in range(max_iter):
# ステップ2: 閾値処理 — 小さい係数をゼロにする
small_inds = np.abs(xi) < threshold
xi[small_inds] = 0
# ステップ3: 残った非ゼロ要素について再度最小二乗
big_inds = ~small_inds
if np.sum(big_inds) == 0:
break
xi[big_inds] = np.linalg.lstsq(
Theta[:, big_inds], dXdt_col, rcond=None
)[0]
return xi
# 各状態変数に対してSTLSを適用
n_vars = X.shape[1]
Xi = np.zeros((Theta.shape[1], n_vars))
threshold = 0.1 # 閾値の設定
for k in range(n_vars):
Xi[:, k] = stls(Theta, dXdt[: , k], threshold=threshold)
# 発見された方程式を表示
var_names = ['dx/dt', 'dy/dt', 'dz/dt']
print("=" * 60)
print("SINDyによって発見された方程式:")
print("=" * 60)
for k in range(n_vars):
equation = f"{var_names[k]} = "
terms = []
for j in range(len(labels)):
if np.abs(Xi[j, k]) > 1e-10:
coeff = Xi[j, k]
terms.append(f"{coeff:+.4f} * {labels[j]}")
equation += " ".join(terms) if terms else "0"
print(equation)
# 真の係数との比較
print("\n" + "=" * 60)
print("真のLorenz方程式の係数との比較:")
print("=" * 60)
print("dx/dt = -10x + 10y")
print("dy/dt = 28x - y - xz")
print("dz/dt = xy - 8/3 z")
# ===== 5. 係数の復元精度の検証 =====
# 真の係数行列を構成
true_Xi = np.zeros_like(Xi)
# dx/dt = -10x + 10y → labels: ['1','x','y','z','xx','xy','xz','yy','yz','zz',...]
x_idx = labels.index('x')
y_idx = labels.index('y')
z_idx = labels.index('z')
xy_idx = labels.index('xy')
xz_idx = labels.index('xz')
true_Xi[x_idx, 0] = -10.0 # dx/dt の x 係数
true_Xi[y_idx, 0] = 10.0 # dx/dt の y 係数
true_Xi[x_idx, 1] = 28.0 # dy/dt の x 係数
true_Xi[y_idx, 1] = -1.0 # dy/dt の y 係数
true_Xi[xz_idx, 1] = -1.0 # dy/dt の xz 係数
true_Xi[xy_idx, 2] = 1.0 # dz/dt の xy 係数
true_Xi[z_idx, 2] = -8.0/3 # dz/dt の z 係数
# 係数の比較
print("係数復元精度の詳細:")
print(f"{'関数':<10} {'変数':<8} {'推定値':<12} {'真値':<12} {'絶対誤差':<12}")
print("-" * 54)
for k in range(n_vars):
for j in range(len(labels)):
if np.abs(Xi[j, k]) > 1e-10 or np.abs(true_Xi[j, k]) > 1e-10:
err = np.abs(Xi[j, k] - true_Xi[j, k])
print(f"{labels[j]:<10} {var_names[k]:<8} {Xi[j,k]:>+10.6f} {true_Xi[j,k]:>+10.6f} {err:>10.2e}")
# 全体の相対誤差
nonzero_mask = np.abs(true_Xi) > 1e-10
relative_error = np.abs(Xi[nonzero_mask] - true_Xi[nonzero_mask]) / np.abs(true_Xi[nonzero_mask])
print(f"\n非ゼロ係数の平均相対誤差: {np.mean(relative_error):.2e}")
print(f"非ゼロ係数の最大相対誤差: {np.max(relative_error):.2e}")
# ===== 6. ノイズを加えた場合の性能検証 =====
noise_levels = [0.0, 0.01, 0.05, 0.1, 0.2]
results = {}
fig, axes = plt.subplots(1, len(noise_levels), figsize=(20, 4))
for idx, noise_level in enumerate(noise_levels):
# ノイズ付加
X_noisy = X + noise_level * np.std(X, axis=0) * np.random.randn(*X.shape)
# 数値微分（ノイズがある場合はSavitzky-Golayフィルタを使うのが理想だが、
# ここでは中心差分で比較）
dXdt_noisy = numerical_derivative(X_noisy, dt)
# ライブラリ構成
Theta_noisy, _ = build_library(X_noisy, poly_order=3)
# STLS適用
Xi_noisy = np.zeros((Theta_noisy.shape[1], n_vars))
for k in range(n_vars):
Xi_noisy[:, k] = stls(Theta_noisy, dXdt_noisy[:, k], threshold=0.1)
# 係数誤差の計算
coeff_error = np.linalg.norm(Xi_noisy - true_Xi, 'fro') / np.linalg.norm(true_Xi, 'fro')
n_nonzero = np.sum(np.abs(Xi_noisy) > 1e-10)
results[noise_level] = {'error': coeff_error, 'n_nonzero': n_nonzero}
# 非ゼロパターンの可視化
axes[idx].imshow(np.abs(Xi_noisy) > 1e-10, aspect='auto', cmap='Blues')
axes[idx].set_title(f'Noise={noise_level*100:.0f}%\nErr={coeff_error:.4f}')
axes[idx].set_xlabel('State variable')
axes[idx].set_ylabel('Library function')
plt.suptitle('SINDy Sparsity Pattern under Different Noise Levels', fontsize=14)
plt.tight_layout()
plt.show()
# 結果のサマリ
print("\nノイズレベルと係数復元精度:")
print(f"{'ノイズ(%)':<12} {'相対誤差':<15} {'非ゼロ項数(真=7)':<20}")
print("-" * 47)
for noise, res in results.items():
print(f"{noise*100:<12.0f} {res['error']:<15.6f} {res['n_nonzero']:<20d}")
# ===== 7. 発見された方程式によるシミュレーション =====
def sindy_simulate(t, state, Xi, labels, var_names_state=['x', 'y', 'z']):
"""SINDyで発見された方程式によるシミュレーション"""
x_val, y_val, z_val = state
X_point = np.array([[x_val, y_val, z_val]])
Theta_point, _ = build_library(X_point, poly_order=3)
dxdt = Theta_point @ Xi
return dxdt.flatten()
# SINDyモデルによるシミュレーション
sol_sindy = solve_ivp(
lambda t, state: sindy_simulate(t, state, Xi, labels),
t_span, x0, t_eval=t_eval, method='RK45',
rtol=1e-10, atol=1e-10
)
# 比較プロット
fig, axes = plt.subplots(3, 1, figsize=(12, 9))
var_labels = ['x(t)', 'y(t)', 'z(t)']
for i in range(3):
axes[i].plot(t, X[:, i], 'b-', label='True', alpha=0.7, linewidth=1)
axes[i].plot(sol_sindy.t, sol_sindy.y[i], 'r--', label='SINDy', alpha=0.7, linewidth=1)
axes[i].set_ylabel(var_labels[i], fontsize=12)
axes[i].legend(fontsize=10)
axes[i].grid(True, alpha=0.3)
axes[0].set_title('Lorenz System: True vs SINDy Reconstruction', fontsize=14)
axes[2].set_xlabel('Time', fontsize=12)
plt.tight_layout()
plt.show()
# 軌道の3D比較
fig = plt.figure(figsize=(12, 5))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot(X[:, 0], X[:, 1], X[:, 2], 'b-', lw=0.5, alpha=0.7)
ax1.set_title('True Lorenz')
ax1.set_xlabel('x'); ax1.set_ylabel('y'); ax1.set_zlabel('z')
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot(sol_sindy.y[0], sol_sindy.y[1], sol_sindy.y[2], 'r-', lw=0.5, alpha=0.7)
ax2.set_title('SINDy Reconstruction')
ax2.set_xlabel('x'); ax2.set_ylabel('y'); ax2.set_zlabel('z')
plt.tight_layout()
plt.show()
# ===== 8. LASSOとの比較 =====
from sklearn.linear_model import Lasso
# LASSO回帰
Xi_lasso = np.zeros_like(Xi)
alpha_lasso = 0.01 # LASSO正則化パラメータ
for k in range(n_vars):
lasso = Lasso(alpha=alpha_lasso, fit_intercept=False, max_iter=10000)
lasso.fit(Theta, dXdt[:, k])
Xi_lasso[:, k] = lasso.coef_
# STLS vs LASSO の係数比較
print("STLS vs LASSO の係数比較:")
print(f"{'関数':<10} {'変数':<8} {'STLS':<12} {'LASSO':<12} {'真値':<12}")
print("-" * 54)
for k in range(n_vars):
for j in range(len(labels)):
if (np.abs(Xi[j,k]) > 1e-10 or np.abs(Xi_lasso[j,k]) > 1e-10
or np.abs(true_Xi[j,k]) > 1e-10):
print(f"{labels[j]:<10} {var_names[k]:<8} "
f"{Xi[j,k]:>+10.6f} {Xi_lasso[j,k]:>+10.6f} {true_Xi[j,k]:>+10.6f}")
# 全体の相対誤差比較
err_stls = np.linalg.norm(Xi - true_Xi, 'fro') / np.linalg.norm(true_Xi, 'fro')
err_lasso = np.linalg.norm(Xi_lasso - true_Xi, 'fro') / np.linalg.norm(true_Xi, 'fro')
print(f"\n全体の相対Frobenius誤差:")
print(f" STLS: {err_stls:.6f}")
print(f" LASSO: {err_lasso:.6f}")
まとめ
本記事では、SINDy（Sparse Identification of Nonlinear Dynamics）について理論から実装まで解説しました。
- SINDyの基本原理: 支配方程式のスパース性を利用し、候補関数ライブラリ $\bm{\Theta}(\bm{X})$ からデータ駆動で方程式を発見する
- STLSアルゴリズム: 反復的な閾値処理と最小二乗法の組み合わせにより、バイアスのないスパース回帰を実現する
- 数値微分の重要性: ノイズに対する感度が高く、適切な微分推定法の選択が精度を左右する
- LASSOとの違い: STLSは係数バイアスがなく、物理法則の正確な係数復元に適している
- PDE-FINDへの拡張: 偏微分方程式の発見にも同じフレームワークが適用可能
- Lorenz方程式の再構成: ノイズがない場合、真の係数をほぼ完全に復元できることを確認した
次のステップとして、以下の記事も参考にしてください。