「線形代数の半歩先」をPythonで写経 ～ ７章 クープマン作用素理論、動的モード分解、未来予測

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

書籍の著者　大久保 潤 先生

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

第５部は「時間発展方程式」を取り扱います。
前回記事から引き続き、 動的モード分解（DMD） の実践です。

第25話～第27話の集大成である「未来予測」をChatGPTに導かれてチャレンジします！

数学素人なので、どうぞお手柔らかにお願いいたします。
では書籍とChatGPTを開いて線形代数の旅に出発です🚀

はじめに

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

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

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

１．前章の振り返り

６章最後尾から振り返ります。

◆ ◆ ◆

そういえば、コード生成時にChatGPTがこんなことを言っていました。

将来の時間へ予測、ですって！？

◆ ◆ ◆

もちろん将来予測を所望します！
焦って依頼文章が途中で切れてしまいましたｗ

憧れの将来予測が遂にこの手に……

◆ ◆ ◆

この記事は、ChatGPTの力を借りて到達できた「未来予測のコード実行」と「コードの解説」に焦点を当てて書きます。

では、この記事のコードで利用するライブラリのインポートからスタートしましょう。

### インポート

# 数値計算
import numpy as np

# 数学演算用の個別関数
from scipy.integrate import solve_ivp   # 常微分方程式の求解

# 動的モード分解(DMD)
from pydmd import DMD
from pydmd.plotter import plot_summary

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

２．ファン・デル・ポール方程式に DMD を適用して未来予測

DMD は動的モード分解（Dynamic Mode Decomposition）のことです。ChatGPT が書いたコードを手がかりにして、テキストの図 25.1 に似せた「時間範囲 $${t = [0, 3]}$$」の可視化と、時刻 $${t=5}$$ までの未来予測の可視化を実装しました。
３つのチャートを描画します。

# ✅ ファン・デル・ポール方程式に DMD を適用する（未来予測つき）

# --- 必要なライブラリのインポート ---
# import numpy as np
# import matplotlib.pyplot as plt
# from scipy.integrate import solve_ivp
# from pydmd import DMD
# from pydmd.plotter import plot_summary

### ファン・デル・ポール方程式による時系列データの生成とDMDの実行

## ファン・デル・ポール方程式の定義
# 非線形の減衰の強さパラメータ（書籍では 𝜖）
mu = 1.0
# ファン・デル・ポール方程式関数の定義：solve_ivpの目的関数
def vdp(t, y):
    x1, x2 = y
    dx1 = x2
    dx2 = mu * x2 * (1 - x1**2) - x1
    return [dx1, dx2]

## 時系列データの生成：ファン・デル・ポル方程式による
# 設定
t_span = (0, 3)    # 時間範囲
delta_t = 0.1      # 時間間隔 Δt
t_eval = np.arange(*t_span, delta_t)
y0 = [0.2, 0.9]    # データの初期値
# ファン・デル・ポール方程式の求解
sol = solve_ivp(vdp, t_span, y0, t_eval=t_eval)
# 解の取得
X = sol.y          # 解：x1,x2
time = sol.t       # 時間間隔

## DMDの実行
dmd = DMD(svd_rank=2)  # svd_rank: 特異値分解のランク数
dmd.fit(X)
plot_summary(dmd)      # DMDの結果を可視化（追加コード）

## 描画処理
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

## モードのプロット
ax1.plot(np.real(dmd.eigs), np.imag(dmd.eigs), 'o')
ax1.set(title=f'DMD 固有値 (複素平面)\n固有値 {dmd.eigs.round(3)}',
        xlabel="実数部", ylabel="虚数部", ylim=(-0.15, 0.15))

## 再構成結果との比較
# x1の元データのデータ点の描画
ax2.plot(time, X[0], 'o', color='tab:blue', label='$x_1$: 数値解')
# x1のDMDによる再構成結果の点線の描画
ax2.plot(time, dmd.reconstructed_data.real[0], '--', color='navy',
         label='$x_1$: DMD再構成')
# x2の元データのデータ点の描画
ax2.plot(time, X[1], 'o', color='tab:orange', label='$x_2$: 数値解')
# x2のDMDによる再構成結果の点線の描画
ax2.plot(time, dmd.reconstructed_data.real[1], '--', color='red',
         label='$x_2$: DMD再構成')
# 修飾
ax2.set(title='ファン・デル・ポール方程式の数値解と\nDMDによる再構成',
        xlabel='時間')
ax2.legend()


### 未来予測

## 設定
# t0からの予測期間の設定
future_period = 5

## DMD パラメータ付与
# t0の時点の値
dmd.dmd_time['t0'] = 0
# 時点の間隔 Δt
dt = time[1] - time[0]
dmd.dmd_time['dt'] = dt
# 予測する時点
future_time = np.arange(3, future_period, delta_t)
future_steps = len(time) + len(future_time)
# 時点数
dmd.dmd_time['tend'] = future_steps

## 予測の実行：DMDを実行して未来分含む再構成データを再取得
# t0からの予測値の取得
future_data = dmd.reconstructed_data.real

## 描画処理：元データ + 未来予測
# 描画領域の設定
fig, ax = plt.subplots(figsize=(10, 5))
# 描画用の時間軸の作成
total_time = dmd.dmd_timesteps * dt
# x1の元データのデータ点の描画
ax.plot(time, X[0], 'o', label='$x_1$: 数値解', color='tab:blue')
# x1のDMDによる再構成結果の点線の描画
ax.plot(total_time, future_data[0], '--', label='$x_1$: DMD予測値', color='navy')
# x2の元データのデータ点の描画
ax.plot(time, X[1], 'o', label='$x_2$: 数値解', color='tab:orange')
# x2のDMDによる再構成結果の点線の描画
ax.plot(total_time, future_data[1], '--', label='$x_2$: DMD予測値', color='red')
# 予測期間最終時点の垂直点線の描画
ax.axvline(future_period, color='gray', linestyle=':')
# 修飾
ax.set(xlabel='時間', ylabel='状態',
       title=f'DMD による未来予測 (t={future_period} まで)')
ax.legend()
ax.grid(True)
plt.show()

【実行結果】
■ DMD のサマリープロット
１行２列目の Discrete-time Eigenvalues で固有値の状況を確認します。
２つの固有値（赤い点）は単位円（点線）付近に存在します。
ChatGPTによると「振幅が保たれる」（減衰しない）を意味するそうです。

■ 時系列データとDMD再構成のプロット
右側のチャートには、ファン・デル・ポール方程式による時系列データが点で示され、時系列データをDMDで近似したデータが点線で示されています。
いい感じの近似に見えます。

■ DMDによる未来予測
遂に時間発展方程式から始まる「未来予測」に到達しました！
予測値を示す２つの点線はいい感じの曲線です！

時間範囲 $${[0, 3]}$$ の時系列データで学習済みのDMDモデル dmd に対して、未来予測のコードで以下のパラメータを与えて取得した「再構築値 reconstructed_data」が予測値に相当するようです。

• 最初の時刻の値 $${t_0=0}$$

• 時刻の間隔$${\Delta t = 0.1}$$

• 未来を含めた計算対象の時点数 引数 tend

たぶん、時刻 $${t=3}$$ 以降の将来を予測するのではなく、DMD が学習した近似関数で時刻 $${t=0}$$ 以降のすべての時点を計算しているようです。
時系列データの生成元であるファン・デル・ポール方程式の近似解を将来時点まで計算している、こんな感じです（あくまで個人の感想です）。

なにはともあれ、第５部の最優先課題「時間発展方程式とそのデータを使って未来予測がしたい！」を成就できました。
やったね！

テキスト第27話最終ページに次のように書かれた「動的モード分解」をPythonで体感できました。

ん？拡張動的モード分解？

ChatGPTの解説 ～教科書～

ここからは、ChatGPTの解説を掲載いたします。
クープマンの作用素理論に基づく動的モード分解（DMD）の初歩的な解説であり、私にとっての「教科書」です！

【おしながき】

1. ファン・デル・ポールの時系列に DMD を適用する（未来予測つき）
   - ChatGPTのオリジナルコードです。

2. DMD 未来予測：コード実行結果の見るべきポイント
   - 上記1. の実行結果に関する解説です。

3. クープマン理論とDMD（基礎から体系的理解へ）
   - 上記1. のコードを以下の観点で解説したものです。

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

１．✅ ファン・デル・ポール方程式に DMD を適用する（未来予測つき）

ChatGPTのオリジナルコードは次の時間を扱っています。

• 時系列データを 時刻 $${t=20}$$ まで生成してDMDで学習

• 未来予測は 時刻 $${t=30}$$ まで

# ✅ ファン・デル・ポール方程式に DMD を適用する（未来予測つき）

--- 必要なライブラリのインポート ---
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from pydmd import DMD

# --- ファン・デル・ポール方程式の定義 ---
mu = 1.0

def vdp(t, y):
    x1, x2 = y
    dx1 = x2
    dx2 = mu * x2 * (1 - x1**2) - x1
    return [dx1, dx2]

# --- 時系列データの生成 ---
t_span = (0, 20)
t_eval = np.linspace(*t_span, 1000)
y0 = [0.2, 0.8]
sol = solve_ivp(vdp, t_span, y0, t_eval=t_eval)

X = sol.y                 # shape = (2, 1000)
time = sol.t              # shape = (1000,)

# --- DMD の適用 ---
dmd = DMD(svd_rank=2)
dmd.fit(X)

# --- 結果の表示 ---
print("DMD 固有値:", dmd.eigs)

# --- モードのプロット ---
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(np.real(dmd.eigs), np.imag(dmd.eigs), 'o')
plt.title("DMD 固有値 (複素平面)")
plt.xlabel("Re")
plt.ylabel("Im")
plt.grid(True)

# --- 再構成結果との比較 ---
plt.subplot(1, 2, 2)
plt.plot(time, X[0], label="$x_1$: original", alpha=0.7)
plt.plot(time, dmd.reconstructed_data.real[0], '--', label="$x_1$: DMD")
plt.plot(time, X[1], label="$x_2$: original", alpha=0.7)
plt.plot(time, dmd.reconstructed_data.real[1], '--', label="$x_2$: DMD")
plt.title("DMD による再構成と元データの比較")
plt.xlabel("Time")
plt.legend()
plt.tight_layout()
plt.show()

# --- 未来予測の設定 ---
future_time = np.linspace(20, 30, 500)
dmd.dmd_time['t0'] = 0
future_steps = len(time) - 1 + len(future_time)
dmd.dmd_time['tend'] = future_steps

dt = time[1] - time[0]
dmd.dmd_time['dt'] = dt

# --- 再構成データの再取得（未来分含む） ---
future_data = dmd.reconstructed_data.real

# --- 時間軸の作成 ---
total_time = np.linspace(0, 30, future_data.shape[1])

# --- 可視化：元データ + 未来予測 ---
plt.figure(figsize=(10, 5))
plt.plot(time, X[0], label="$x_1$: original", color='C0')
plt.plot(time, X[1], label="$x_2$: original", color='C1')
plt.plot(total_time, future_data[0], '--', label="$x_1$: DMD prediction", 
         color='C0')
plt.plot(total_time, future_data[1], '--', label="$x_2$: DMD prediction",
         color='C1')
plt.axvline(20, color='gray', linestyle=':')
plt.title("DMD による未来予測 (t=20 以降)")
plt.xlabel("Time")
plt.ylabel("State")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

【実行結果】

時系列データ（実線）は時間 $${t=3}$$ までのチャートと比べて、波形が乱れている印象です。
一方で予測値（点線）は時間経過とともに上下の振幅が大きくなっています。

続いて、このチャートの結果に関するChatGPTの解説に移りましょう。

２．🧠 DMD 未来予測：コード実行結果の見るべきポイント

🔍 出力グラフの構成
DMD による未来予測のグラフは以下のような構成になっています：

• 横軸：時間（例：0 〜 30秒）

• 縦軸：状態変数 $${x_1}$$, $${x_2}$$

• 実線：元データ（0〜20秒）

• 破線：DMD による再構成＋未来予測（0〜30秒）

• 点線の縦線：予測の開始点（$${t = 20}$$）

✅ チェックポイント

$$
\begin{array}{l:l}
項目 & 観察点と解釈 \\
\hline
 \\
元データ \text{vs DMD}再構成 & 一致具合、周期、振幅 \\
& →\text{DMD} モデルがどれだけ忠実に\\
&　元データを再現しているか \\
 \\
未来予測の滑らかさ & 崩壊・発散していないか \\
& →\text{DMD}モデルの安定性を示す \\
 \\
周期の継続性 & 周期が保たれているか \\
& →固有値の虚部が非ゼロで \\
& 　振動モードが捉えられている証拠 \\
 \\
振幅の維持 & 急激に大きく/小さくなっていないか \\
& →固有値の絶対値 ≈ 1\ は \\
& 　増幅も減衰もしないことを意味 \\
\end{array}
$$

🧮 固有値の解釈（例）

たとえば DMD 固有値が：

$$
\lambda = 1.00067 \pm 0.01987i
$$

このとき：

• $${|\lambda| \approx 1}$$ → 振幅が保たれる

• $${\text{Im}(\log(\lambda)) \neq 0}$$ → 振動成分あり

周期は：

$$
T = \frac{2\pi}{\text{Im}(\log(\lambda))} \approx \frac{2\pi}{0.01987} \approx 316
$$

※ 単位は時間ステップ幅（`dt`）に依存

✅ 理想的な出力はこうなる！

$$
\begin{array}{l:l}
条件 & 理由 \\
\hline
 \\
再構成（破線）が & モデル構築が成功している \\
元データとよく一致 & \\
 \\
未来予測（t>20）も & クープマンモードが物理的に妥当 \\
滑らかで周期が継続 & \\
 \\
振幅が発散せず安定 & 固有値の絶対値 ≈ 1 を保っている \\
\end{array}
$$

🔁 さらに進めるなら

• 固有値の時間スケール変換（周期・減衰定数などの数値計算）

• クープマンモードの空間的な構造可視化

• ランクを変えての再実験と比較分析

• 将来予測精度を評価する指標（MSE等）

DMD は時系列データに対する強力なツールなので、どんどん実験してみましょう！

３．✅ クープマン理論とDMD（基礎から体系的理解へ）

未来予測コードを動かした後、テキストの主要な用語とコードとの関わりを知りたくなって、尋ねました。

以下がChatGPTの回答＝教科書です。

◆ ◆ ◆

🧠 クープマン理論とDMD：概念マップ

コルモゴロフ後退方程式（連続時間の観測量の期待値のPDE）
       ↓（離散化）
クープマン作用素 K（無限次元の線形演算子）
       ↓（有限次元に近似）
クープマン行列 A ≈ K（観測量の行列による線形近似）
       ↓（固有分解）
固有値 λ、固有関数 φ、モード（ベクトルの方向成分）
       ↓（時間発展の再構成）
g(x(t)) ≈ Σ モード × e^{λt}

◆ ◆ ◆

🔍 各概念の役割と使われ方

$$
\begin{array}{l:l}
概念 & コードやモデル内での役割 \\
\hline
 \\
①コルモゴロフ & 観測値の期待値を時間的に予測するための \\
　後退方程式 & 偏微分方程式（\text{PDE}） \\
&→ \text{DMD}で使うことは少ないが、 \\
&　連続時間のクープマン作用素の \\
&　微分形として重要 \\
 \\
②クープマン作用素 & 観測関数の時間発展を表す線形作用素 \\
&→ \text{DMD}はこの作用素を有限次元で \\
&　近似しようとしている \\
 \\
③クープマン行列 & クープマン作用素の有限次元近似。 \\
 （\text{Koopman Matrix}） & \text{DMD}ではこれを実際に計算している  \\
 \\
④クープマンモード & 観測量の空間的な特徴を表す成分。 \\
& \text{DMD}で得られるモードベクトル \\
 \\
⑤クープマン固有値 & クープマン行列の固有値 \\
&（時間方向の周期・固有関数・減衰）、 \\
& 固有ベクトル・関数 \\
&（変化しない特徴）\\
 \\
⑥クロネッカー積 & 複数変数の高次特徴量の拡張 \\
& （辞書関数）やテンソル化に使われる \\
& （例：\text{EDMD}など）
\end{array}
$$

◆ ◆ ◆

① コルモゴロフ後退方程式：連続時間版の基礎

数式：

$$
\cfrac{\partial u}{\partial t} = \mathcal{L} u
$$

• $${u(x, t) = \mathbb{E}[g(X_t) \mid X_0 = x]}$$

• $${\mathcal{L}}$$ は拡散過程の生成演算子

関係：

• DMDはこの PDE を 「離散時間で近似的に解く」 アプローチと解釈できる

• よって、DMD ≒ 後退方程式の有限次元近似

◆ ◆ ◆

② クープマン作用素：観測量の時間発展を司る本体

数式定義：

$$
(Kg)(x)=\mathbb{E}[g(f(x))]またはg(x(t+1))=Kg(x(t))
$$

関係：

• DMDはこの作用を 行列 $${A}$$ によって近似

• 実装例

X1 = X[:, :-1]                            # g(x_t)
X2 = X[:, 1:]                             # g(x_{t+1})
# A ≈ X2 @ X1.T @ (X1 @ X1.T)^-1

• この $${A}$$ がクープマン作用素 $${K}$$ の有限次元近似（＝クープマン行列）

◆ ◆ ◆

③ クープマン行列：クープマン作用素の具体的な行列表現

• DMDでは、クープマン作用素 $${K}$$ の有限次元近似を行います。
  それが：

$$
A \approx K
$$

• この $${A}$$ が、DMDの中で クープマン行列（Koopman matrix）と呼ばれるものです
  実際にDMDで求めているのがこの「近似行列」$${A}$$です

• Eigen-decomposition に使われる対象です

• pydmd.DMD では、クープマン行列自体は直接アクセスできる属性はありません

• 実装例：
  このnumpyのSVD実装で得られる A_tilde が、DMDで求められる クープマン行列の有限次元近似です

# X1: 過去の観測系列（時刻 t〜t+T-1）
# X2: 次の観測系列（時刻 t+1〜t+T）

X1 = X[:, :-1]
X2 = X[:, 1:]

# SVD 分解
U, S, Vh = np.linalg.svd(X1, full_matrices=False)
S_inv = np.diag(1 / S)

# クープマン行列の近似
A_tilde = X2 @ Vh.T @ S_inv @ U.T

◆ ◆ ◆

④ クープマンモード：空間方向の“かたち”

• 固有ベクトルに対応する観測空間の方向ベクトル

• 観測値の空間的構造を抽出する

• 実装例：

    modes = dmd.modes  # shape: (n_features, n_modes)

• 描画や解釈対象となる

◆ ◆ ◆

⑤ クープマン固有値・固有関数

• 固有値 $${\lambda_i}$$：時間変化のスケール（周期、減衰）

• 固有関数：$${K\phi = \lambda \phi}$$ を満たす観測量（非線形でも時間とともに形が変わらない）

• 実装例：

    eigvals = dmd.eigs

• 時間発展を $${e^{\lambda t}}$$ でモデル化している部分！

◆ ◆ ◆

⑥ クロネッカー積：高次特徴量の拡張（EDMDなど）

• 状態が $${(x_1, x_2)}$$ だとすると、辞書関数 $${\phi(x) = [x_1, x_2, x_1^2, x_1 x_2, x_2^2]}$$ のような構造を作る必要がある

• そのときに、クロネッカー積が自然な拡張として使われる

• 実装例：

import numpy as np
x = np.array([1, 2])
phi = np.kron(x, x)  # → [1, 2, 2, 4] = [x1*x1, x1*x2, x2*x1, x2*x2]

• EDMD（Extended DMD）で明示的に使われるが、DMDの拡張として理解すると役立つ

◆ ◆ ◆

✅ まとめ：コード実装との対応

$$
\begin{array}{l:l}
概念 & 実装での役割、関係の深い箇所 \\
\hline
 \\
コルモゴロフ後退方程式 & 背景理論（連続時間） \\
& \text{PDE}に基づく意味付け \\
 \\
クープマン作用素 & 時間発展の本体 \\
& X[:,1:] = A @ X[:,:-1] \\
 \\
クープマン行列 & \text{DMD}の目的 \\
& \text{A\_tilde}（近似行列）\\
 \\
クープマンモード & 空間構造の可視化 \\
& \text{dmd.modes} \\
 \\
固有値・関数 & 周期・減衰の指標 \\
& \text{dmd.eigs} \\
 \\
クロネッカー積 & 特徴拡張・非線形性対応 \\
& \text{EDMD}, カーネル\text{DMD}などで使用 \\
\end{array}
$$

おわりに

「まとめ」の内容が気になります。

ChatGPTは常に興味を惹く提案をしてくれます。
で、聞いてみました。

この回答にはさまざまなテーマが連なっていて、どれをやっていいのか判断がつきませんでした。
そこで…

DMD 完全マスターの道ですって！

次章へ続く！

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

１．のんびり統計

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

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

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

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

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

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

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

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

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

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

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

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

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

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