カルマンスムーザ3種比較（RTS固定区間・固定ラグ・固定点）：理論・遅延と計算量のトレードオフ・Python実装

はじめに：フィルタリングとスムージングの違い

状態空間モデルにおける状態推定には、使える観測の範囲によって 2 つの問題設定があります。

• フィルタリング：時刻 \(k\) までの観測 \(y_{1:k}\) を用いた推定 → \(p(x_k \mid y_{1:k})\) 。因果的（causal） な処理で、観測が届くたびに推定を更新できるため オンライン・リアルタイム処理 に適します。
• スムージング：時刻 \(k\) より未来の観測も用いた推定 → \(p(x_k \mid y_{1:T})\) （\(T > k\) ）。非因果的（non-causal） な処理で、未来の情報を使うぶんフィルタより常に同等以上の精度が得られますが、その分の 遅延 が発生します。

「未来の観測をどこまで待つか」によって、スムーザは大きく 3 種類に分類されます。

1. 固定区間スムーザ（fixed-interval smoother）：全観測 \(y_{1:T}\) が揃ってから、全時刻 \(k = 0, \ldots, T\) の \(\hat{x}_{k|T}\) を推定する。代表が Rauch-Tung-Striebel（RTS）スムーザ。完全オフライン。
2. 固定ラグスムーザ（fixed-lag smoother）：常に \(N\) ステップ前の状態 \(\hat{x}_{k-N|k}\) を推定する。遅延 \(N\) を許容する準リアルタイム処理。
3. 固定点スムーザ（fixed-point smoother）：特定の時刻 \(j\) の状態 \(\hat{x}_{j|k}\) を、観測が増えるたびに更新し続ける。初期状態やイベント発生時刻の推定に使う。

本記事では 3 種のスムーザの漸化式を導出し、計算量・遅延・用途を比較したうえで、numpy フルスクラッチの Python 実装でカルマンフィルタ・固定ラグ・RTS の RMSE を実測比較します。

基盤となる理論は カルマンフィルタの理論とPython実装 を参照してください。

3種スムーザの比較表

はじめに結論の比較表を示します。\(n\) は状態次元、\(T\) はデータ長、\(N\) はラグ長です。

カルマンフィルタの復習（前向きパス）

3 種のスムーザはいずれもカルマンフィルタの前向きパスの結果を入力とします。線形ガウス状態空間モデルを定義します。

予測ステップ

更新ステップ

スムージングでは各時刻の \(\hat{x}_{k|k}, P_{k|k}, \hat{x}_{k|k-1}, P_{k|k-1}\) を保存しておく必要があります。非線形系では式 \((3)\) 〜\((7)\) を EKF や UKF の対応するステップに置き換えることで、同じ枠組みのままスムーザを構成できます。

固定区間スムーザ（RTS）

固定区間スムージングは、全観測 \(y_{1:T}\) を使って全時刻の事後分布 \(p(x_k \mid y_{1:T})\) を求める問題です。RTS スムーザは時刻 \(T\) から \(0\) へ向かう後向き再帰で解きます。

スムーザゲイン：

平滑化状態・共分散：

初期条件は \(\hat{x}_{T|T}, P_{T|T}\) （前向きパスの最終結果）で、\(k = T-1, \ldots, 0\) の順に適用します。\(P_{k+1|T} \preceq P_{k+1|k}\) より \(P_{k|T} \preceq P_{k|k}\) が常に成り立ち、スムーザの不確実性はフィルタ以下になります。導出の詳細は RTSスムーザの記事 を参照してください。

固定ラグスムーザ

固定ラグスムーザは、時刻 \(k\) の観測を受け取るたびに \(N\) ステップ前の平滑化推定 \(\hat{x}_{k-N|k}\) を出力します。ストリーミング処理のまま「\(N\) ステップ分だけ未来を見た」推定が得られるのが特徴です。

実現方法1：拡大状態カルマンフィルタ

過去 \(N+1\) 時刻分の状態を縦に積んだ拡大状態ベクトル

を定義し、対応する拡大遷移行列（1 ブロックずつシフトする構造）に対して通常のカルマンフィルタをそのまま回す方法です（Biswas–Mahalanabis 固定ラグスムーザ）。実装は単純ですが、状態次元が \((N+1)n\) になるため計算量は \(O((Nn)^3)\) とラグ長に対して急速に増えます。

実現方法2：窓付きRTS（本記事の実装）

各時刻 \(k\) で、直近の窓 \([k-N, k]\) に限定して RTS の後向き再帰（式 \((8)\) 〜\((9)\) ）を回し、窓の先頭 \(\hat{x}_{k-N|k}\) だけを出力する方法です。1 ステップあたりの計算量は \(O(N n^3)\) で済み、前向きパスの中間結果を直近 \(N\) 個だけ保持すればよいためメモリも \(O(N n^2)\) です。

ラグ長 \(N\) と精度のトレードオフ

未来の観測がフィルタ推定に与える修正は、スムーザゲイン \(G_k\) を通じて後向きに伝搬します。定常状態でのゲインを \(\bar{G}\) とすると、\(N\) ステップ先の観測の寄与はおおむね \(\bar{G}^N\) のオーダーで減衰します。つまり：

• \(N\) を増やすほど RTS（\(N = \infty\) 相当）の精度に漸近するが、改善は指数的に飽和する
• 遅延とメモリは \(N\) に線形に増加する

実務では「\(\bar{G}^N\) が十分小さくなる \(N\) 」（システムの時定数の数倍程度）を選べば、RTS とほぼ同等の精度が遅延 \(N\) で得られます。後述の実験でこの飽和挙動を数値で確認します。

固定点スムーザ

固定点スムーザは、関心のある特定の時刻 \(j\) （例：初期時刻、イベント発生時刻）の推定 \(\hat{x}_{j|k}\) を、新しい観測 \(y_k\) （\(k > j\) ）が届くたびに更新します。

時刻 \(j\) から \(k-1\) までのスムーザゲインの積を

と定義すると、フィルタのイノベーション \(\nu_k = y_k - H \hat{x}_{k|k-1}\) とカルマンゲイン \(K_k\) を用いて

という 1 ステップ \(O(n^3)\) の逐次更新が得られます。時刻 \(j\) の情報だけを追跡するためメモリは \(O(n^2)\) で固定です。\(G_i\) の積 \(B_k\) は \(k\) が進むにつれ縮小していくため、更新量は次第に小さくなり、\(\hat{x}_{j|k}\) は固定区間スムーザの \(\hat{x}_{j|T}\) に漸近します。

Python実装：ランダムウォークでのRMSE比較

1 次元ランダムウォーク＋観測ノイズという最もシンプルな系で、カルマンフィルタ・固定ラグスムーザ・RTS スムーザの精度を比較します。

モデル定義とデータ生成

import numpy as np

# ---- 状態空間モデル：ランダムウォーク + 観測ノイズ ----
T = 300  # タイムステップ数

F = np.array([[1.0]])   # 状態遷移（ランダムウォーク）
H = np.array([[1.0]])   # 観測行列
Q = np.array([[0.05]])  # プロセスノイズ分散
R = np.array([[1.0]])   # 観測ノイズ分散

n, m = 1, 1
x0 = np.array([0.0])
P0 = np.diag([1.0])

カルマンフィルタ（前向きパス）

def kalman_filter(y_obs, F, H, Q, R, x0, P0):
    """カルマンフィルタの前向きパス（スムーザ用に中間結果をすべて保存）"""
    T = len(y_obs)
    n = len(x0)
    x_filt = np.zeros((T + 1, n))       # x_{k|k}
    P_filt = np.zeros((T + 1, n, n))    # P_{k|k}
    x_pred = np.zeros((T, n))           # x_{k|k-1}
    P_pred = np.zeros((T, n, n))        # P_{k|k-1}
    x_filt[0], P_filt[0] = x0, P0
    for k in range(T):
        # 予測（式3, 4）
        x_pred[k] = F @ x_filt[k]
        P_pred[k] = F @ P_filt[k] @ F.T + Q
        # 更新（式5, 6, 7）
        S = H @ P_pred[k] @ H.T + R
        K = P_pred[k] @ H.T @ np.linalg.inv(S)
        x_filt[k + 1] = x_pred[k] + K @ (y_obs[k] - H @ x_pred[k])
        P_filt[k + 1] = (np.eye(n) - K @ H) @ P_pred[k]
    return x_filt, P_filt, x_pred, P_pred

RTSスムーザ（後向きパス）

def rts_smoother(x_filt, P_filt, x_pred, P_pred, F):
    """RTS（固定区間）スムーザの後向きパス"""
    T = len(x_pred)
    n = x_filt.shape[1]
    x_smooth = np.zeros((T + 1, n))
    P_smooth = np.zeros((T + 1, n, n))
    x_smooth[T], P_smooth[T] = x_filt[T], P_filt[T]
    for k in range(T - 1, -1, -1):
        G = P_filt[k] @ F.T @ np.linalg.inv(P_pred[k])  # 式8
        x_smooth[k] = x_filt[k] + G @ (x_smooth[k + 1] - x_pred[k])  # 式9
        P_smooth[k] = P_filt[k] + G @ (P_smooth[k + 1] - P_pred[k]) @ G.T  # 式10
    return x_smooth, P_smooth

固定ラグスムーザ（窓付きRTS）

def fixed_lag_smoother(y_obs, F, H, Q, R, x0, P0, N):
    """固定ラグスムーザ：各時刻 k で直近 N ステップ分だけ後向き再帰を回し
    x_{k-N|k} を出力する（窓付きRTSによる実装）"""
    T = len(y_obs)
    n = len(x0)
    x_filt, P_filt, x_pred, P_pred = kalman_filter(y_obs, F, H, Q, R, x0, P0)
    x_lag = np.zeros((T + 1, n))
    x_lag[0] = x0
    for k in range(1, T + 1):
        j0 = max(0, k - N)  # 出力する時刻 k-N
        xs = x_filt[k].copy()
        for j in range(k - 1, j0 - 1, -1):  # 窓 [j0, k] で後向き再帰
            G = P_filt[j] @ F.T @ np.linalg.inv(P_pred[j])
            xs = x_filt[j] + G @ (xs - x_pred[j])
        x_lag[j0] = xs
    # 末尾 N 点はラグ分の未来観測が無いので、最終窓の後向き再帰で埋める
    xs = x_filt[T].copy()
    x_lag[T] = xs
    for j in range(T - 1, max(0, T - N) - 1, -1):
        G = P_filt[j] @ F.T @ np.linalg.inv(P_pred[j])
        xs = x_filt[j] + G @ (xs - x_pred[j])
        x_lag[j] = xs
    return x_lag

モンテカルロ評価（50試行）

単一系列だとサンプリングの揺らぎで順位が入れ替わることがあるため、numpy.random.default_rng で 50 系列を生成して平均 RMSE を比較します。

def rmse(a, b):
    return np.sqrt(np.mean((a - b) ** 2))

n_trials = 50
lags = [1, 2, 3, 5, 10, 20]
rmse_kf_list, rmse_rts_list = [], []
rmse_lag_lists = {N: [] for N in lags}

rng = np.random.default_rng(42)
for trial in range(n_trials):
    # 真の状態と観測の生成
    x_true = np.zeros((T + 1, n))
    y_obs = np.zeros((T, m))
    for k in range(T):
        x_true[k + 1] = F @ x_true[k] + rng.multivariate_normal(np.zeros(n), Q)
        y_obs[k] = H @ x_true[k + 1] + rng.multivariate_normal(np.zeros(m), R)

    x_filt, P_filt, x_pred, P_pred = kalman_filter(y_obs, F, H, Q, R, x0, P0)
    x_rts, P_rts = rts_smoother(x_filt, P_filt, x_pred, P_pred, F)

    rmse_kf_list.append(rmse(x_true[1:, 0], x_filt[1:, 0]))
    rmse_rts_list.append(rmse(x_true[1:, 0], x_rts[1:, 0]))
    for N in lags:
        x_lag = fixed_lag_smoother(y_obs, F, H, Q, R, x0, P0, N)
        rmse_lag_lists[N].append(rmse(x_true[1:, 0], x_lag[1:, 0]))

rmse_kf = np.mean(rmse_kf_list)
rmse_rts = np.mean(rmse_rts_list)
print(f"KF  filter RMSE : {rmse_kf:.4f}")
for N in lags:
    r = np.mean(rmse_lag_lists[N])
    print(f"Fixed-lag N={N:2d}   : {r:.4f}  ({(1 - r / rmse_kf) * 100:.1f}% improvement)")
print(f"RTS smoother    : {rmse_rts:.4f}  ({(1 - rmse_rts / rmse_kf) * 100:.1f}% improvement)")

実行結果

KF  filter RMSE : 0.4519
Fixed-lag N= 1   : 0.4145  (8.3% improvement)
Fixed-lag N= 2   : 0.3890  (13.9% improvement)
Fixed-lag N= 3   : 0.3717  (17.7% improvement)
Fixed-lag N= 5   : 0.3521  (22.1% improvement)
Fixed-lag N=10   : 0.3397  (24.8% improvement)
Fixed-lag N=20   : 0.3387  (25.0% improvement)
RTS smoother    : 0.3388  (25.0% improvement)

期待どおり フィルタ < 固定ラグ < RTS の順で精度が改善し、固定ラグは \(N\) の増加とともに RTS に漸近します（\(N=20\) で差はモンテカルロ誤差の範囲内）。

この系の定常予測分散は \(\bar{P}_{k|k-1} = 0.25\) で、定常スムーザゲインは \(\bar{G} = \bar{P}_{k|k} / \bar{P}_{k+1|k} = 0.8\) です。未来観測の寄与は \(\bar{G}^N = 0.8^N\) で減衰するため、\(N=10\) で \(0.8^{10} \approx 0.11\) 、\(N=20\) で \(\approx 0.012\) となり、実測の「\(N=10\) でほぼ飽和、\(N=20\) で RTS と同等」という結果と整合します。遅延をどこまで許容できるかでラグ長を決め、\(\bar{G}^N\) が数%を切る \(N\) を上限の目安にするのが実務的な設計指針です。

発展：非線形・非ガウスへの拡張

本記事は線形ガウス系に限定しましたが、同じ 3 分類は非線形・非ガウス系にもそのまま持ち上がります。EKF/UKF ベースのスムーザに加え、粒子フィルタ の重み付きサンプルを後向きに再重み付けする粒子スムーザ（forward-filtering backward-smoothing、\(O(T M^2)\) 、\(M\) は粒子数）を使えば、多峰性のある事後分布の平滑化も可能です。さらに近年は、状態遷移・観測モデル自体をニューラルネットで学習する Deep State Space Model（Deep Kalman Filter、S4/Mamba 系列など）が発展しており、変分推論で近似した事後分布に対して同様のフィルタリング／スムージングの区別が現れます。古典スムーザの遅延・計算量の設計感覚は、これらの現代的手法を使う際にもそのまま役立ちます。

まとめ

• スムージングは未来の観測を使う非因果的な推定で、フィルタリングより常に同等以上の精度が得られる代わりに遅延が生じる
• RTS（固定区間） は全時刻を最高精度で推定するオフラインの標準解、固定ラグ は遅延 \(N\) と引き換えに準リアルタイムでスムージング精度を得る折衷案、固定点 は特定時刻だけを逐次精密化する軽量な方式
• ランダムウォーク系の 50 試行モンテカルロで、フィルタ RMSE 0.4519 → 固定ラグ \(N=5\) で 0.3521 → RTS で 0.3388 と単調に改善することを確認した
• ラグ長 \(N\) の効果は定常スムーザゲイン \(\bar{G}\) の \(N\) 乗で指数的に飽和するため、システムの時定数の数倍を目安に選べばよい

おすすめ書籍

カルマンフィルタの導出から実装例までを日本語で丁寧に解説した定番書です。状態空間モデルの理解が深まります。

※ 上記は Amazon アソシエイトのリンクです。

関連記事

• カルマンスムーザ（RTS Smoother）の理論・前向き後向き再帰・Python実装 - 固定区間スムーザ（RTS）の導出と実装を詳説しています。本記事の前提となる記事です。
• カルマンフィルタの理論とPython実装 - 全スムーザの前向きパスとなるカルマンフィルタの基礎です。
• 拡張カルマンフィルタ（EKF）の理論とPython実装 - 非線形系への第一歩。EKF ベースの非線形スムーザの基盤です。
• Unscented Kalman Filter（UKF）の理論とPython実装 - シグマポイントによる非線形推定。UKF スムーザへ拡張できます。
• 粒子フィルタのPython実装：リサンプリング手法の比較 - 非ガウス系の粒子スムーザの基盤となる粒子フィルタを解説しています。

参考文献

• Rauch, H. E., Tung, F., & Striebel, C. T. (1965). “Maximum likelihood estimates of linear dynamic systems.” AIAA Journal, 3(8), 1445-1450.
• Biswas, K. K., & Mahalanabis, A. K. (1972). “An approach to fixed-lag smoothing problems.” IEEE Transactions on Aerospace and Electronic Systems, AES-8(5), 676-682.
• Meditch, J. S. (1969). Stochastic Optimal Linear Estimation and Control. McGraw-Hill.
• Särkkä, S. (2013). Bayesian Filtering and Smoothing. Cambridge University Press.