「線形代数の半歩先」をPythonで写経～９章コルモゴロフ後退方程式、DMD時系列クラスタリング

ネイピア DS

2025年4月21日 11:51

第５部「ならべた数のさらなる発展」

書籍の著者　大久保潤先生

この記事は、書籍「線形代数の半歩先」の第５部「ならべた数のさらなる発展」に掲載の「非線形系における線形性」に関する Python写経活動のドキュメンタリーです。

第５部は「時間発展方程式」を取り扱います。

前回に引き続き、ChatGPT提案の「完全マスター」を実践します。
ChatGPTによるコードと解説に全面依拠して 動的モード分解（DMD）マスターを目指してまいります。

コルモゴロフ後退方程式 と DMD特徴量を活用した時系列クラスタリング を実践します！

いっけんテキストから離れるように見えますが、動的モード分解の理解を通じて、テキストを読み解く素地ができると思います！
とにかく数学素人なので、どうぞお手柔らかにお願いいたします。

では書籍とChatGPTを開いて線形代数の旅に出発です🚀

はじめに

このブログシリーズは、書籍「線形代数の半歩先データサイエンス・機械学習に挑む前の30話」（講談社サイエンティフィク、「テキスト」と呼びます）の Python 写経の実践を通じて得た個人的な知見を書きます。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

第５部ならべた数のさらなる発展

クープマン理論＆DMD完全マスター・ロードマップ

ChatGPT提案の「クープマン理論＆DMD完全マスター」ロードマップです。（DMD は動的モード分解の略称です）

この記事は STEP 4 と STEP 5 を実践します。
主に、コルモゴロフ後退方程式（オルンシュタイン・ウーレンベック過程）、DMD特徴量を活用した時系列クラスタリングに取り組みます。

テキストとの大凡の関連を表にしました。

$$
\begin{array}{l:l}
テーマ & テキスト関連箇所 \\
\hline
\\
コルモゴロフ後退方程式 & \text{p.235} \\
\end{array}
$$

（注意事項）
ChatGPTの回答をそのまま記載しています。
内容の適否はチェックしていませんので、ご了承ください。

STEP4: 生成演算子の数値近似（コルモゴロフ後退方程式の離散化）

オルンシュタイン・ウーレンベック過程（OU過程）に対して、コルモゴロフ後退方程式の数値解法を用いる（L行列による差分法）

（１）Python 実装
コルモゴロフ後退方程式の数値解法の実装です。
将来時刻 $${T}$$ に区間 $${[-1, 1]}$$ にいる確率を可視化します。

# --- セル 1: ライブラリとグリッド定義 ---
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # ★はマニュアル追加コード 
from scipy.sparse import diags
from scipy.sparse.linalg import expm

# グリッドとパラメータ
x_min, x_max = -5, 5
N = 500
x = np.linspace(x_min, x_max, N)
dx = x[1] - x[0]

# OU過程パラメータ（Kolmogorovの後退方程式形式）
gamma = 1.0  # ドリフト（復元力）係数
D = 0.5      # 拡散係数（= σ² / 2）

# --- セル 2: 安定な差分法による生成演算子Lの定義（中心差分・Neumann境界） ---
# ドリフト項：-gamma * x * ∂u/∂x
c = -gamma * x
main1 = np.zeros(N)
off1 = c[1:] / (2 * dx)
L_drift = diags([-off1, main1, off1], offsets=[-1, 0, 1], shape=(N, N)).toarray()

# 拡散項：D * ∂²u/∂x²（Neumann境界：端の2点は固定）
main2 = -2 * D / dx**2 * np.ones(N)
off2 = D / dx**2 * np.ones(N-1)
L_diff = diags([off2, main2, off2], offsets=[-1, 0, 1], shape=(N, N)).toarray()

# 総合演算子 L
L = L_drift + L_diff

# --- セル 3: 境界の調整（Neumann条件を明示的に） ---
L[0] = np.zeros(N)
L[-1] = np.zeros(N)
L[0, 0] = -1.0  # u(x_min) = u(x_min)
L[-1, -1] = -1.0

# --- セル 4: 観測関数（最終時間での初期値） ---
indicator = np.where((x >= -1) & (x <= 1), 1.0, 0.0)

# --- セル 5: 後退時間発展 ---
T_back = 1.0
u_xt = expm(T_back * L) @ indicator

# --- セル 6: 可視化 ---
plt.figure(figsize=(8, 5))
plt.plot(x, u_xt, label=r"$u(x, T)$")
plt.title("コルモゴロフ後退方程式による予測（OU過程 / 安定差分法）")
plt.xlabel("x（現在の位置）")
plt.ylabel("確率（将来に [-1,1] にいる確率）")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

【実行結果】
両裾が切り立っている感じが面白いですね。

■ Python コードに書かれていた概要説明

🔧 ポイント

差分法：中心差分（一次）＋ Neumann 境界（端点固定）
行列Lを非対称・発散しにくく調整
行列指数も SciPy sparse版（`scipy.sparse.linalg.expm`）使用で高速化

✅ 今後の応用先

固有関数の可視化（Lの固有分解）
ファン・デル・ポールへの拡張
後退方程式とクープマン固有分解の接続

（２）🔍 グラフの読み解き

📉 横軸：現在の位置 $${x}$$
📈 縦軸：将来時刻 $${T}$$ において区間 $${[-1, 1]}$$ にいる確率

🔍 解釈のポイント

$$
\begin{array}{l:l}
観察 & 意味 \\
\hline
\\
左右対称な滑らかな曲線 & OU過程はガウス的で \\
& 対称な復元力（ドリフト）を持つ \\
\\
x=0 のとき最大 & 中心（平均）にいるほど、\\
& 将来 [−1,1] に留まる確率が最大 \\
\\
遠ざかるにつれ & 現在地が ±5 など遠いと、そこから \\
なだらかに減少 & [−1,1] に到達する確率が低くなる \\
\\
確率値は 0〜1 の範囲内に & 発散や負の値なし \\
収まっている & → 数値的にも安定 \\
\end{array}
$$

🧠 知識として押さえておきたい背景

OU過程の解析的な確率密度はガウス分布
後退方程式での $${u(x,T)}$$ は「未来条件に従った関数値の期待値」
指示関数を使ったので、それは「将来ある区間にいる確率」

🎯 総合評価：この出力は OK？

$$
\begin{array}{l:l:l}
項目 & 評価 & コメント \\
\hline
\\
対称性 & ✅ & 物理的に妥当 \\
範囲の収まり & ✅ & 0〜1 の確率として正しい \\
数値安定性 & ✅ & 発散なし、なめらか \\
解釈しやすさ & ✅ & 確率として直感的にわかりやすい \\
\end{array}
$$

正しく実装されており、物理・数理的に非常に自然です！

（３）✅ コルモゴロフ後退方程式（OU過程）の数値解法から学べる注目ポイント・要点整理

1️⃣ コルモゴロフの後退方程式とは？

🔍 数式の構造
後退方程式の一般形：

$$
\cfrac{\partial u}{\partial t} =Lu
$$

OU過程における生成演算子：

$$
Lu (x) = - \gamma x \cfrac{\partial u}{\partial x} + D \cfrac{\partial^2u}{\partial x^2} \\
$$

💡 学び

時間の逆方向に進む偏微分方程式
$${u(x,T)}$$ ：将来における観測値（例：確率）を現在位置 $${x}$$ から見たときの予測量
期待値計算の枠組みとして、マルコフ過程と深く結びついている

2️⃣ 数値的な生成演算子 $${L}$$ の構築

🔍 実装ポイント

# ドリフト項（一次微分）：中心差分
# 拡散項（二次微分）：Neumann境界（反射）

💡 学び

偏微分項（∂/∂x, ∂²/∂x²）を差分法で行列に近似
安定性を保つには、中心差分やNeumann境界が重要
$${L}$$ は非対称な疎行列になるが、固有分解や指数関数で扱える

3️⃣ 観測関数の考え方：未来条件を初期値に使う

🔍 実装ポイント

indicator = np.where((x >= -1) & (x <= 1), 1.0, 0.0)

💡 学び

$${u(x,T)}$$ は未来における「ある事象が起こるか」の指示関数
それを現在に向かって後退させることで、任意の $${x}$$ に対して「その未来が起こる確率」が得られる

4️⃣ 行列指数関数による時間発展

🔍 実装ポイント

u_xt = expm(L * T_back) @ indicator

💡 学び

時間発展の本質は $${u(t) = e^{Lt} u(T)}$$
行列の指数関数（exponential）で「連続時間の解」が得られる
時間離散化せずに一気に未来（または過去）にジャンプできる

5️⃣ グラフの解釈：現在の位置と将来確率の関係

🔍 実装ポイント

plt.plot(x, u_xt)

💡 学び

グラフは「今この場所にいたら、将来指定した事象が起こる確率」を表している
これは予測関数であり、生成演算子 $${L}$$ の支配構造を可視化したもの
対称性・滑らかさ・最大位置などから物理的/確率的直感が得られる

🎓 このノートブックで得られるスキルセット

$$
\begin{array}{l:l}
分野 & 内容 \\
\hline
\\
数理モデル & \text{OU}過程・後退方程式の構造理解 \\
数値解析 & 差分法・行列指数・境界条件処理 \\
プログラミング & \text{NumPy / SciPy / Sparse}行列の扱い \\
可視化と解釈 & 確率予測関数の意味を読み解く力 \\
\end{array}
$$

🔜 応用への展開アイデア

$$
\begin{array}{l:l}
次にやること & 学べること \\
\hline
\\
✅ L の固有分解と可視化 & クープマン固有関数との関係 \\
✅ 他の観測関数（例：x^2, \sin x）& 関数の進化と性質の比較 \\
✅ ファン・デル・ポール方程式 & 非線形系での後退方程式の \\
& 構造への拡張 \\
\end{array}
$$

STEP5: DMDモードに基づく時系列クラスタリング（t-SNE/PCA可視化 + k-means）

DMDを活用した時系列データ分析に関しては、ChatGPTと試行錯誤しながら学びを進めました。
その結果、クラスタリングをはじめとする「４つのテーマ」がまとまりました。
テーマと見出しの関係を整理します。

【STEP 5 の目次】

１．時系列クラスタリング
（１）Pyhon 実装
（２）グラフの読み解き
（３）注目ポイント・要点整理

２．近傍時系列を抽出して可視化
（４）Python 実装

３．DMD で生み出す時系列データ特徴量の解説
（５）DMDによる「時系列の特徴量化」完全ガイド

４．DMDモードの振幅による寄与の違いの可視化
（６）Python実装
（７）注目ポイント・要点整理

（１）Python 実装
時系列データのクラスタリングを実装します。
30 個の時系列データに DMD を適用して、t-SNE と k-means 法でクラスタリングします。

# --- セル 1: ライブラリとデータ準備（複数時系列を合成） ---
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # ★はマニュアル追加コード 
from pydmd import DMD
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# シンプルな3種の系（正弦波 + 減衰系 + 非線形風）
def gen_sin(freq, t):
    return np.sin(2 * np.pi * freq * t)

def gen_decay(freq, t):
    return np.exp(-0.2 * t) * np.cos(2 * np.pi * freq * t)

def gen_combo(freq, t):
    return np.sin(2*np.pi*freq*t) * (1 - np.exp(-0.5*t))

# 時系列データseriesの生成
T = np.linspace(0, 10, 200)
series = []
labels_true = []

for i in range(10):
    # 正弦波, ラベル:0
    series.append(gen_sin(1 + 0.2*i, T))
    labels_true.append(0)
    # 減衰系, ラベル:1
    series.append(gen_decay(1 + 0.2*i, T))
    labels_true.append(1)
    # 非線形風, ラベル:2
    series.append(gen_combo(1 + 0.2*i, T))
    labels_true.append(2)

series = np.array(series)            # shape (30, 200)　30系列、200時点
labels_true = np.array(labels_true)  # 正解ラベル

# --- セル 2: 各時系列に対して DMD を適用し、モード特徴を抽出 ---
dmd_features = []  # 特徴量：固有値（実数部）, 固有値（虚数部）, 振幅

for s in series:
    dmd = DMD(svd_rank=5)
    dmd.fit(s.reshape(1, -1))
    eigvals = dmd.eigs[:5]
    amps = np.abs(dmd.amplitudes[:5])
    features = np.hstack([np.real(eigvals), np.imag(eigvals), amps])
    dmd_features.append(features)

dmd_features = np.array(dmd_features) 

# --- セル 3: 特徴量の標準化と t-SNE 次元削減 ---
# 標準化
scaler = StandardScaler()
features_scaled = scaler.fit_transform(dmd_features)
# t-SNEによる次元削減の実行（2次元に圧縮） ★ perplexityを20から21へ変更
tsne_proj = TSNE(n_components=2, perplexity=21, random_state=42
            ).fit_transform(features_scaled)

plt.figure(figsize=(6, 5))
plt.scatter(tsne_proj[:, 0], tsne_proj[:, 1], c=labels_true, cmap='Set1', s=80,
            alpha=0.5)
plt.title("DMDモードに基づく時系列の埋め込み可視化（t-SNE）")
plt.xlabel("Dim 1")
plt.ylabel("Dim 2")
plt.grid(True)
plt.tight_layout()
plt.show()

# --- セル 4: k-means によるクラスタリング ---
kmeans = KMeans(n_clusters=3, random_state=0)
kmens_preds = kmeans.fit_predict(features_scaled)

plt.figure(figsize=(6, 5))
plt.scatter(tsne_proj[:, 0], tsne_proj[:, 1], c=kmens_preds, cmap='Set1', s=80,
            alpha=0.5)
plt.title("DMDクラスタリング結果（k-means）")
plt.xlabel("Dim 1")
plt.ylabel("Dim 2")
plt.grid(True)
plt.tight_layout()
plt.show()

【実行結果】

DMD の「固有値・実部」「固有値・虚部」「振幅」を特徴量とし、t-SNE で２次元に次元圧縮して、散布図をプロットしています。
点の色は「正解ラベル」です。
赤とグレイが混じっている点がいくつか見られます。
30 個の時系列データに対して、散布図にあらわれているのは 24 個。
6 個の点が重なっています。

こちらは、３つの特徴量を用いて、k-means 法で３つのクラスタに分けて、t-SNE で圧縮した２次元の散布図にしています。
点の色は k-means 法によるクラスタを示します。
重なっている点は同じクラスタに分けられた模様です。

時系列データといえば、横軸に時間を置いて、時系列の推移を見たいものです。
追加コードをChatGPTに作ってもらったので、後でご紹介します。

■ Python コードに書かれていた概要説明

🔧 ポイント

特徴量に DMD固有値（実部＋虚部）に加えてモード振幅（amplitudes）も追加
特徴ベクトルを `StandardScaler` で正規化（t-SNEやk-meansの性能向上）
t-SNEの `perplexity=20` に設定し、近傍関係の精度向上
※注：perplexity=21 に変更して実行しています。

✅ 効果

クラスタ分離性の改善
可視化の2次元性の復活（横一直線の回避）
DMDの動的構造に基づくクラスタリングの実用化

（２）🔍 グラフの読み解き

✅ 上段（t-SNE可視化）について

$$
\begin{array}{l:l}
観察 & 解釈 \\
\hline
\\
明確に2つのグループに & \text{DMD}の特徴ベクトルが \\
分かれている & 系の違いを反映している証拠 \\
\\
\text{t-SNE}が2次元できれいな & スケーリング + \text{perplexity} \\
分布を描いている & 調整が機能している \\
\\
点は滑らかに分布しており、 & 入力特徴量が正しく \\
極端な外れ値がない & 圧縮されている \\
\end{array}
$$

💡 意味

t-SNE空間での分布は「DMDモード特徴量の距離構造」を反映しています。
時系列が、潜在的な力学的類似性で埋め込まれていると解釈できます。

✅ 下段（k-meansクラスタリング）について

$$
\begin{array}{l:l}
観察 & 解釈 \\
\hline
\\
クラスタが3つに & \text{n\_clusters=3}の指定通り \\
分かれている & 機能している \\
\\
一部の領域（中央）で色の & 境界的な例がグレーゾーンで \\
混ざりがある & 分類されている可能性 \\
\\
左右の分布は上段の & \text{k-means}のクラスタリング結果が \\
可視化結果と整合している & 視覚的にも理解しやすい状態にある \\
\end{array}
$$

💡 意味

k-meansクラスタリングは、DMD固有値＋モードのパターンに基づいて時系列を分類できている。
この手法は、ラベルなし時系列データのクラスタリングや異常検知にも応用可能です。

🏁 総合評価：この出力は妥当か？

$$
\begin{array}{l:l:l}
評価軸 & 判定 & コメント \\
\hline
\\
特徴抽出 & ✅ & \text{DMD}によって動的情報が \\
& & 抽出できている \\
\\
次元削減 & ✅ & \text{t-SNE}が有効に機能し、 \\
& & 構造が埋め込まれている \\
\\
クラスタリング & ✅ & \text{k-means}で分離されており、 \\
& & 探索的分析に適している \\
\\
可視化の直感性 & ✅ & 分布と分類結果が視覚的に \\
& & わかりやすい \\
\end{array}
$$

→ 非常に妥当で、成功したクラスタリング結果です！

（３）✅ DMDモードに基づく時系列クラスタリング：注目ポイント・要点整理

1️⃣ DMD は時系列の“力学的な特徴抽出器”

💡 要点

各時系列に対して DMD を適用すると：

固有値（eigenvalues）：時間周波数・成長減衰を表す
振幅（amplitudes）：そのモードがどれだけ支配的か
時系列を構成するダイナミクスを、少数の数値（スペクトル）で表現できる

これが「時系列の特徴量化」に最適

2️⃣ DMD固有値・モードを「特徴量ベクトル」に変換

features = [Re(λ1), ..., Im(λ1), ..., |amp1|, ...]

💡 要点

複素固有値 → 実部（成長/減衰）＋虚部（周波数成分）
モード振幅 → ダイナミクスの重要度
これらを連結して、特徴ベクトル（DMDスペクトル）として使用

3️⃣ 特徴空間の次元削減で「似た時系列」を可視化

💡 要点

SNE（または PCA）を使って 2次元に埋め込むことで、

時系列の類似性を可視化
潜在構造（クラスタ）を発見できる
DMD特徴は物理的・力学的に意味のある座標なので、時系列の“力学的空間”上で比較できる

4️⃣ クラスタリング（k-means）による分類

💡 要点

DMD特徴に基づいてクラスタリングすることで：

“振る舞いの似ている時系列”をグループ化できる
ラベルなしの時系列でも意味のある分類が可能に
時系列分類・異常検知・プロセス監視に有効

5️⃣ 正常に機能させるための工夫ポイント

$$
\begin{array}{l:l}
工夫 & 意図・効果 \\
\hline
\\
✅ 特徴量にモード振幅を追加 & 固有値だけでは足りない情報を補う \\
✅ \text{StandardScaler}で & \text{t-SNE}や\text{k-means}に必要な前処理 \\
スケーリング & \\
✅ \text{perplexity} の調整 & \text{t-SNE} の分離性能を向上 \\
\end{array}
$$

🎓 このノートブックで得られるスキルセット

$$
\begin{array}{l:l}
分野 & 内容 \\
\hline
\\
時系列解析 & \text{DMD}を使った力学的表現の獲得 \\
\\
機械学習 & 特徴抽出 → 次元削減 \\
& → クラスタリングのパイプライン構築 \\
\\
可視化 & \text{t-SNE/PCA}での潜在空間表現 \\
\\
実践力 & ノイズ・振動・医療データ等への \\
& 応用の発想力 \\
\end{array}
$$

🔜 応用・発展のアイデア

$$
\begin{array}{l:l}
アイデア & 内容 \\
\hline
\\
📈 各クラスタの代表時系列を比較 & 系の特徴の視覚的理解 \\
\\
📊 クラスタごとの & 力学構造の違いを解析 \\
\text{DMD} 固有値分布を比較 & \\
\\
✅ 外れ値検知 & 特徴空間で孤立する点を \\
& 異常検知として使う \\
\\
🔁 動画・振動・センサーデータ & 汎用的な異常検知・ \\
への応用 & 行動分類へ展開可能\\
\end{array}
$$

✨ まとめ

DMDは「時系列に潜むダイナミクス」を特徴量として取り出し、クラスタリングを通じて「力学的に似た系列」を見つけ出せる強力な方法です！

（４）📈 （追加コード）近傍時系列を抽出して可視化

30 個の時系列データから「似ている」ものを抽出して、時系列推移をプロットします。
２個目の時系列データ（インデックスは $${1}$$）と３つのDMD特徴量「固有値・実部」「固有値・虚部」「振幅」の「ユークリッド距離」が近い５つの時系列データを描画します。

# --- セル 1〜4 までそのまま（前回のクラスタリングまで） ---
# ※ 省略：ここまでで proj, features_scaled, series, preds, labels_true 
# が得られている前提

# --- セル 5: 近傍系列を取得する関数 ---
from sklearn.metrics.pairwise import euclidean_distances

def find_similar_series(index, features_scaled, preds, proj, topk=5):
    dist_feat = euclidean_distances(features_scaled[index:index+1],
                                    features_scaled).flatten()
    dist_tsne = euclidean_distances(proj[index:index+1], proj).flatten()

    similar_same_cluster = np.where(preds == preds[index])[0]
    nearest_by_distance = np.argsort(dist_feat)[1:topk+1]  # 自身は除外
    nearest_by_tsne = np.argsort(dist_tsne)[1:topk+1]

    return {
        '同クラスタ内': similar_same_cluster.tolist(),
        '特徴空間の近傍': nearest_by_distance.tolist(),
        't-SNE空間の近傍': nearest_by_tsne.tolist(),
    }

# --- セル 6: 対象インデックスと近傍を取得 ---
target_index = 1  # 比較したい時系列のインデックス
neighbors = find_similar_series(
    target_index, features_scaled, kmens_preds, tsne_proj, topk=5)

# --- セル 7: 可視化 ---
plt.figure(figsize=(10, 6))
plt.plot(T, series[target_index], label=f"Target ({target_index})", lw=3,
         color='black')

for i, idx in enumerate(neighbors['特徴空間の近傍']):
    plt.plot(T, series[idx], label=f"Neighbor {i+1} (idx={idx})", lw=0.8)

plt.title("近傍時系列の比較（DMD特徴空間ベース）")
plt.xlabel("Time")
plt.ylabel("Signal")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

【実行結果】
６個の時系列データは、同時刻では異なる位置にいますが、減衰の傾向がとても似ています。

■ find_similar_series() 関数の徹底解説
このコードの核心とも言える、近似時系列データを特定する関数「find_similar_series()」について解説してもらいました。

✅ 関数の目的

ある時系列の「似ている系列」を複数の観点から探し出すための関数です。
3種類の類似度ベースで「近い時系列のインデックスリスト」を返します：

$$
\begin{array}{l:l}
類似度の種類 & 意味 \\
\hline
\\
クラスタ\text{ID}ベース & 同じクラスタに属する系列 \\
& （粗い分類）\\
\text{DMD}特徴空間ベース & \text{DMD}スペクトルで近い\\
& （固有値＋モード振幅） \\
\text{t-SNE}空間ベース & 2次元埋め込み上で近い\\
& （直感的距離）
\end{array}
$$

🔍 関数の定義

def find_similar_series(index, features_scaled, preds, proj, topk=5):

$$
\begin{array}{l:l}
引数 & 意味 \\
\hline
\\
\text{index} & 比較したいターゲット時系列のインデックス \\
\text{features\_scaled} & \text{DMD}から得た特徴量（標準化済み） \\
\text{preds} & クラスタリング結果（\text{k-means}の出力） \\
\text{proj} & \text{t-SNE}などの2次元埋め込み結果 \\
\text{topk} & 近い系列を何件取得するか（デフォルトは5件）\\
\end{array}
$$

🔧 内部処理

① DMD特徴空間での距離を計算

dist_feat = euclidean_distances(features_scaled[index:index+1],
                                features_scaled).flatten()

features_scaled は各時系列の DMD特徴ベクトル
index 番目の特徴量ベクトルと他のすべてのベクトルとのユークリッド距離を計算

② t-SNE空間での距離を計算

dist_tsne = euclidean_distances(proj[index:index+1], proj).flatten()

2次元に圧縮された空間（t-SNE座標）で、同様に距離を測る

③ 同クラスタ内の系列を抽出

similar_same_cluster = np.where(preds == preds[index])[0]

ターゲット系列と同じクラスタIDを持つ系列すべてを抽出

④ 近い系列を距離順に取得（自身を除く）

nearest_by_distance = np.argsort(dist_feat)[1:topk+1]
nearest_by_tsne = np.argsort(dist_tsne)[1:topk+1]

np.argsort() で距離が小さい順に並べ、最も近い上位 topk 件を抽出
[1:] にしているのは、自分自身（距離0）を除くため

✅ 戻り値（辞書）

return {
    '同クラスタ内': similar_same_cluster.tolist(),
    '特徴空間の近傍': nearest_by_distance.tolist(),
    't-SNE空間の近傍': nearest_by_tsne.tolist(),
}

🔙 返ってくるのは、3種類のインデックスリスト

$$
\begin{array}{l:l}
キー & 内容 \\
\hline
\\
同クラスタ内 & \text{preds} で同じクラスタ\text{ID}を持つ全系列 \\
特徴空間の近傍 & \text{DMD}特徴ベクトルで近い順に \text{topk} 件 \\
\text{t-SNE}空間の近傍 & 埋め込み空間上で近い順に \text{topk} 件 \\
\end{array}
$$

💡 この関数でできること

似た系列を見つける手段を切り替えて比較できる
- DMDスペクトルで「力学的に」似ているか
- t-SNEで「視覚的に」近いか
- クラスタIDで「ざっくり分類的に」近いか
異常点やクラスタ境界の分析にも使える
- 類似点が全然見つからない点 → 外れ値の候補

🎯 実践応用の例

$$
\begin{array}{l:l}
使い方 & 応用シーン \\
\hline
\\
ユーザの行動履歴の類似系列推薦 & レコメンド \\
センサーデータの異常検出 & 故障予知・診断 \\
複数センサの時系列構造比較 & 系のクラスタ解析 \\
\end{array}
$$

（５）✅ DMDによる「時系列の特徴量化」完全ガイド

🔍 そもそも「特徴量化」とは？

生の時系列（例：200次元）を少数の意味ある数値ベクトル（例：10次元）に変換すること
通常の方法：
- 統計特徴（平均・分散など）
- 周波数特徴（FFT、スペクトログラム）
- 埋め込み特徴（Autoencoderなど）

📌 DMDでの特徴量化は？

時系列の時間変化を「振動＋増減」の要素に分解し、その要素を特徴量にすること！

🧠 概念ステップ

時系列を以下のようにみなします：

$$
x(t) \approx \sum_{k=1}^r a_k \phi_k e^{\omega_k t}
$$

$$
\begin{array}{l:l}
要素 & 意味 \\
\hline
\\
\phi_k & 固有モード（空間的形） \\
\omega_k & 固有値（振動・成長率）\\
a_k & 振幅（寄与度）\\
\end{array}
$$

これらの数値的情報（固有値・振幅など）をベクトル化することで、時系列を特徴量化できます。

🛠 実装ステップ

✅ 1. 時系列を行列化（DMD形式に整形）

x = s.reshape(1, -1)  # shape = (1, T) にする（PyDMDに合わせる）

✅ 2. DMDを適用

from pydmd import DMD

dmd = DMD(svd_rank=5)
dmd.fit(x)

svd_rank=5：上位5モードだけ抽出（次元削減にもなる）

✅ 3. 固有値を抽出して特徴に使う

eigvals = dmd.eigs[:5]  # 複素数（λ₁, λ₂, ..., λ₅）

複素数をそのまま使うのではなく、以下のように分解して使います：

real_parts = np.real(eigvals)
imag_parts = np.imag(eigvals)

✅ 4. 振幅（モードの寄与度）を加える

amps = np.abs(dmd.amplitudes[:5])

振幅が大きいほど、そのモードが時系列の変化に貢献している

✅ 5. 特徴ベクトルにまとめる

feature_vec = np.hstack([real_parts, imag_parts, amps])

結果：1つの時系列が 15次元のベクトルに変換される（5モード × 3成分）

🎯 得られる特徴量ベクトルの意味

$$
\begin{array}{l:l}
要素 & 意味 \\
\hline
\\
実部（Re(λ)） & 成長率・減衰率（正→発散、負→収束） \\
虚部（Im(λ)） & 振動数（周期性） \\
振幅（amplitude） & モードの寄与度 \\
\end{array}
$$

→ これらを組み合わせて、「その時系列がどんな動的性質を持っているか」を数値で表現できます！

✅ 振幅（amplitude, |a|）の意味

DMDで得られる各モードには、そのモードがどれだけ時系列の構造に寄与しているかを示す係数「振幅（amplitude）」が対応します。

この振幅は、初期条件（時系列の最初の状態）を、DMDモードの線形結合で再構成する際の係数に相当します。
そのため、振幅が大きいモードほど、その時系列に強く現れている重要な動的パターンであると解釈できます。

振幅が大きい → そのモードが時系列の波形や変動に対して大きな影響を与えている
振幅が小さい → 時系列にほとんど現れていないモード（ノイズ、微弱な成分）

つまり、振幅は「モードの重要度・支配力」そのものです。

💡 使いどころ・応用例

$$
\begin{array}{l:l}
応用 & 内容 \\
\hline
\\
✅ クラスタリング & \text{DMD}特徴ベクトルに基づいて \\
& 動的構造で分類 \\
\\
✅ 検索・推薦 & 似た\text{DMD}特徴の系列を検索して提示 \\
\\
✅ 異常検知 & 特徴空間上で孤立している時系列を \\
& 異常と判定 \\
\\
✅ 予測モデルの説明 & 学習済モデルの特徴空間と \\
&重ね合わせて説明可能に \\
\end{array}
$$

📦 コードまとめ（1つの時系列を特徴ベクトルにする）

def dmd_feature_vector(s, svd_rank=5):
    from pydmd import DMD
    s = s.reshape(1, -1)
    dmd = DMD(svd_rank=svd_rank)
    dmd.fit(s)
    eigvals = dmd.eigs[:svd_rank]
    amps = np.abs(dmd.amplitudes[:svd_rank])
    return np.hstack([np.real(eigvals), np.imag(eigvals), amps])

（６）✅ DMDモードの振幅による寄与の違いを可視化する

# ✅ DMDが複数モードを抽出できるようにする：スナップショット行列の構築と解析

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # ★はマニュアル追加コード 
from pydmd import DMD

# --- 時系列の生成（複雑な信号） ---
T_full = np.linspace(0, 10, 500)

def gen_complex_signal(t):
    return (
        np.sin(2 * np.pi * 1.0 * t) +
        0.5 * np.sin(2 * np.pi * 3.0 * t) +
        0.8 * np.sin(2 * np.pi * 5.0 * t) * np.exp(-0.1 * t)
    )

x_full = gen_complex_signal(T_full)

# --- スナップショット行列の作成（windowed embedding） ---
# 各スナップショットの長さ
window_size = 50
# スナップショットの数
num_snapshots = len(x_full) - window_size
# スナップショット行列Xの作成 shape: (window_size, num_snapshots)
X = np.column_stack([x_full[i:i+window_size] for i in range(num_snapshots)])

# --- DMDの実行 ---
dmd = DMD(svd_rank=10)
dmd.fit(X)
# 固有値eigvalsとモードampsの取得
eigvals = dmd.eigs
amps = np.abs(dmd.amplitudes)
sorted_indices = np.argsort(amps)[::-1]

# --- モード数と特異値を確認 ---
print(f"Number of DMD modes: {dmd.modes.shape[1]}")
print(f"Eigenvalues: \n{eigvals}")

# --- 可視化（複数モードの再構成） ---
t0 = T_full[-window_size:]  # 表示用時間軸（最後のスナップショットに対応）
dt = (T_full[-1] - T_full[0]) / (len(T_full) - 1)

plt.figure(figsize=(10, 6))
plt.plot(t0, X[:, -1], label="Original Snapshot", color='black', lw=2)

for i in range(min(5, len(sorted_indices))):
    idx = sorted_indices[i]
    phi = dmd.modes[:, idx]  # shape = (50,)
    omega = np.log(eigvals[idx]) / dt
    time_dynamics = dmd.amplitudes[idx] * np.exp(omega * (t0 - t0[0])) # shape=(50,)
    recon_i = (phi * time_dynamics).real  # element-wise product
    plt.plot(t0, recon_i, label=f"Mode {i+1} (|a|={amps[idx]:.2f})",
             linestyle='--')

plt.title("各DMDモードの寄与（スナップショット行列使用）")
plt.xlabel("Time")
plt.ylabel("Signal")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

【実行結果】

# DMD サマリーの可視化
from pydmd.plotter import plot_summary
plot_summary(dmd, index_modes=(0, 2, 4),)

【実行結果】

■ Python コードに書かれていた概要説明

🎯 スナップショット行列構築のポイント

単一時系列をウィンドウに分割することで、多様な構造を含む「多次元データ」としてDMDに渡す
これにより、DMDが時間的な変化の構造（振動・減衰）を複数モードとして抽出可能に
特に信号の周波数構成や局所的な減衰・成長の理解に有効

（７）🔍DMDモードの振幅による寄与の違いを可視化するの評価ポイント

✅ モードの波形

5つの異なるDMDモードが、時間的に異なる振動成分をしっかり表現しています。
各モードの周波数・振幅・位相が異なるのが視覚的にわかります。
Mode 1, 2 のように振幅が大きいモードは全体形状に大きく寄与しており、
Mode 4, 5 のような細かい振動は高周波成分の補正を担っているように見えます。

✅ 元信号との一致

黒線（Original Snapshot）とモードの和の形状が非常に近く、
→ DMDによる高精度な分解・近似ができている証拠です。

✅ ラベルや色分け

|a|（振幅）表示もあるため、「どのモードがどれだけ支配的か」も一目瞭然です。
視認性もよく、モード解析の教材や報告資料にそのまま使えるレベルです 💯

💡 次のおすすめアクション

$$
\begin{array}{l:l}
やってみたいこと & 内容 \\
\hline
\\
🔁 累積モード再構成 & \text{Mode 1 → Mode 1+2 → ...} と順に \\
& 足していき「再構成の改善」を確認 \\
\\
📈 周波数 \text{vs} 振幅プロット & 固有値の虚部から周波数を算出し、 \\
& スペクトル的に分析 \\
\\
⏩ 未来予測 & 固有値で時間発展して、未来の信号を \\
& 合成・予測（\text{DMD}の強み） \\
\\
📦 クラスタリングや & 複数時系列で\text{DMD}モードの \\
特徴抽出 & 空間にマッピングして分類や分析 \\
\end{array}
$$

おわりに

DMD による「時間発展的な時系列データ」の分析を満喫できましたね！

そして完全マスターの道は…まだまだ続きます！

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

https://note.com/e_dao/n/na614e7d62d81

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルはＲ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

https://note.com/e_dao/n/n60fe8e5413e1

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

https://note.com/e_dao/n/n8c6c82709541

４．楽しい写経ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

https://note.com/e_dao/n/n44dd7c7a645f

５．RとStanではじめる心理学のための時系列分析入門を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

https://note.com/e_dao/n/n96c820e3cac2

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

https://note.com/e_dao/n/ne28cc4931c79

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

https://note.com/e_dao/n/n782bc4b7c5ba

最後までお読みいただきまして、ありがとうございました。

いいなと思ったら応援しよう！

応援ありがとうございます。これからもがんばって記事を作成します！

この記事が参加している募集

新生活をたのしく

8.3万件

「線形代数の半歩先」をPythonで写経 ～ ９章 コルモゴロフ後退方程式、DMD時系列クラスタリング

第５部「ならべた数のさらなる発展」

はじめに

第５部 ならべた数のさらなる発展

クープマン理論＆DMD完全マスター・ロードマップ

STEP4: 生成演算子の数値近似（コルモゴロフ後退方程式の離散化）

STEP5: DMDモードに基づく時系列クラスタリング（t-SNE/PCA可視化 + k-means）

おわりに

シリーズの記事

ブログの紹介

いいなと思ったら応援しよう！

この記事が参加している募集

「線形代数の半歩先」をPythonで写経～９章コルモゴロフ後退方程式、DMD時系列クラスタリング

第５部ならべた数のさらなる発展