カルマンスムーザ（Rauch-Tung-Striebel/RTS Smoother）の理論・前向き後向き再帰・Python実装（numpy.linalg.inv）

はじめに：フィルタとスムーザの違い

時系列データから隠れた状態を推定する問題を考えます。カルマンフィルタは逐次的にデータを処理し、各時刻で「現在までの観測」を用いて状態を推定します。一方、スムーザは「すべての観測データ」を用いて各時刻の状態を推定します。

• フィルタリング：時刻 \(t\) までの観測 \(y_{1:t}\) を用いた推定 → \(p(x_t | y_{1:t})\)
• スムージング：全時刻の観測 \(y_{1:T}\) を用いた推定 → \(p(x_t | y_{1:T})\)

スムーザは未来の情報も活用できるため、フィルタよりも常に同等以上の推定精度を達成します。ただし、すべてのデータが揃ってから処理するため、リアルタイム処理には適用できません。バッチ処理や事後解析（軌道の再推定、パラメータ推定の前処理など）で威力を発揮します。

本記事では、最も代表的な固定区間スムーザである Rauch-Tung-Striebel（RTS）スムーザ の理論を導出し、Pythonで実装します。

関連記事： 指数移動平均（EMA）フィルタの周波数特性 、 Cubature Kalman Filter（CKF）の理論とPython実装 、 信号処理におけるフィルタリング手法の基礎

カルマンフィルタの復習

RTSスムーザの入力はカルマンフィルタのフォワードパスの出力です。ここで線形ガウス状態空間モデルを定義し、カルマンフィルタの式を簡単に振り返ります。

状態空間モデル

ここで \(F\) は状態遷移行列、\(H\) は観測行列、\(Q\) はプロセスノイズの共分散、\(R\) は観測ノイズの共分散です。

予測ステップ

更新ステップ

RTSスムーザを実行するためには、フォワードパスで各時刻の \(\hat{x}_{k|k}\) 、\(P_{k|k}\) 、\(\hat{x}_{k|k-1}\) 、\(P_{k|k-1}\) をすべて保存しておく必要があります。

RTSスムーザの導出

RTSスムーザは、カルマンフィルタのフォワードパスの結果を用いて、時刻 \(T\) から \(0\) に向かって逆方向に再帰的に平滑化を行います。ここでは「なぜこの後ろ向き再帰が正しい平滑化分布 \(p(x_k \mid z_{1:T})\) を与えるのか」を、ベイズの定理とマルコフ性だけから天下りなしに証明します。

なぜ後ろ向き再帰で正しい平滑化分布が求まるのか

状態空間モデルの同時分布は、状態のマルコフ性（\(x_k\) は \(x_{k-1}\) にのみ依存）と観測の条件付き独立性（\(z_k\) は \(x_k\) にのみ依存）から、次のように分解されます。

平滑化分布 \(p(x_k \mid z_{1:T})\) （\(k<T\) ）を、1 ステップ先の平滑化分布 \(p(x_{k+1} \mid z_{1:T})\) から再帰的に求めるために、まず次の補題を示します。

補題：\(p(x_k \mid x_{k+1}, z_{1:T}) = p(x_k \mid x_{k+1}, z_{1:k})\) （未来の観測 \(z_{k+1:T}\) を知っていても、\(x_{k+1}\) が与えられれば \(x_k\) の条件付き分布は過去の観測 \(z_{1:k}\) だけで決まる）

証明：観測列を \(z_{1:T} = \{z_{1:k}, z_{k+1:T}\}\) と分割し、ベイズの定理を適用すると

（比例定数は \(x_k\) に依存しない正規化項）。マルコフ性より、時刻 \(k+1\) 以降の状態・観測はすべて \(x_{k+1}\) を起点として生成されるため、\(z_{k+1:T}\) は \(x_{k+1}\) が与えられれば \(x_k\) や \(z_{1:k}\) と条件付き独立です。同様に状態遷移も \(x_{k+1}\) は \(x_k\) のみに依存します。

これらを代入すると

先頭の因子 \(p(z_{k+1:T} \mid x_{k+1})\) は \(x_k\) に依存しないため、\(x_k\) についての比例関係を変えません。よって

（最後の式変形は、同じベイズの定理を \(z_{1:T}\) の代わりに \(z_{1:k}\) について適用したものに他なりません）。両辺とも \(x_k\) についての確率密度として正規化されるため、\(p(x_k \mid x_{k+1}, z_{1:T}) = p(x_k \mid x_{k+1}, z_{1:k})\) が成り立ちます。∎

この補題は「未来の観測が \(x_k\) の推定に与える影響は、すべて \(x_{k+1}\) を経由してのみ伝わる」というスムーザの本質を表しています。この補題と周辺化・ベイズの定理から

が得られます。右辺の \(p(x_k \mid x_{k+1}, z_{1:k})\) はフォワードパスの量だけから計算できる分布（次節で導出）、\(p(x_{k+1} \mid z_{1:T})\) は時刻 \(k+1\) まで（\(k+1\) から \(T\) へ向かって）既に計算済みの平滑化分布です。つまりこの式こそが「時刻 \(T\) から \(0\) へ向かう後ろ向き再帰」の数学的根拠です。

平滑化ゲイン \(G_k\) の導出（ガウス分布の条件付け）

\(p(x_k \mid x_{k+1}, z_{1:k})\) を具体的に求めます。線形ガウスモデルでは、時刻 \(k\) までの観測 \(z_{1:k}\) が与えられたとき、フィルタリング分布は \(x_k \mid z_{1:k} \sim \mathcal{N}(\hat{x}_{k|k}, P_{k|k})\) です。また状態遷移式 \((1)\) より \(x_{k+1} = F x_k + w_k\) で、プロセスノイズ \(w_k\) は \(z_{1:k}\) が与えられたときの \(x_k\) と独立です。したがって \((x_k, x_{k+1})\) は \(z_{1:k}\) のもとで同時ガウス分布に従い、平均・共分散は

となります（交差共分散 \(\mathrm{Cov}(x_k, x_{k+1} \mid z_{1:k}) = \mathrm{Cov}(x_k, F x_k + w_k \mid z_{1:k}) = P_{k|k} F^\top\) 、対角ブロックは予測ステップの式 \((3)(4)\) そのものです）。

多変量ガウス分布の条件付け公式（シューア補元）を用いると、同時分布 \(\begin{bmatrix} u \\ v \end{bmatrix} \sim \mathcal{N}\!\left(\begin{bmatrix}\mu_u\\ \mu_v\end{bmatrix}, \begin{bmatrix}\Sigma_{uu} & \Sigma_{uv} \\ \Sigma_{vu} & \Sigma_{vv}\end{bmatrix}\right)\) に対して

が成り立ちます。\(u = x_k\) 、\(v = x_{k+1}\) として代入すると

これがスムーザゲイン \(G_k\) です。\(P_{k|k}F^\top\) は「時刻 \(k\) から見た \(k+1\) への交差共分散」、\(P_{k+1|k}^{-1}\) はそれを予測共分散のスケールで正規化する項と解釈できます。

平滑化状態推定と共分散の導出（周辺化）

上で得た \(p(x_k \mid x_{k+1}, z_{1:k})\) は \(x_{k+1}\) の線形関数を平均に持つガウス分布です。\(p(x_{k+1}\mid z_{1:T}) = \mathcal{N}(\hat{x}_{k+1|T}, P_{k+1|T})\) （時刻 \(k+1\) で既に計算済みの平滑化分布）を用いて前節の積分（周辺化）を実行すると、ガウス線形変換の性質（平均は線形写像をそのまま適用し、共分散には \(x_{k+1}\) 側の不確実性の伝播分が加わる）から

が得られます。これが冒頭で提示したRTSスムーザの式 \((8)\) 〜\((10)\) そのものであり、ベイズの定理・マルコフ性・ガウス分布の条件付けという3つの初等的な事実だけから厳密に導出されたことになります。直感的には、式 \((9)\) は「未来の情報によってフィルタの予測 \(\hat{x}_{k+1|k}\) がどれだけ修正されたか」を、スムーザゲイン \(G_k\) を介して時刻 \(k\) に伝搬しています。

\(P_{k|T} \leq P_{k|k}\) の証明（後ろ向き帰納法）

\(k=T\) から \(k=0\) への逆方向の数学的帰納法で示します。

• base：\(k=T\) では \(P_{T|T} = P_{T|T}\) （平滑化分布とフィルタリング分布が一致し、自明に等号成立）。
• 帰納段階：\(P_{k+1|T} \preceq P_{k+1|k}\) （半正定値順序、すなわち \(P_{k+1|k}-P_{k+1|T}\) が半正定値行列）と仮定すると、式 \((10)\) より

右辺は半正定値行列 \(P_{k+1|k}-P_{k+1|T}\) に対する合同変換 \(G_k(\cdot)G_k^\top\) であり、合同変換は半正定値性を保つため、\(P_{k|k}-P_{k|T}\) も半正定値、すなわち \(P_{k|T} \preceq P_{k|k}\) です。

帰納法により、すべての \(k=0,\ldots,T\) で \(P_{k|T}\preceq P_{k|k}\) 、つまりスムーザの共分散は常にフィルタの共分散以下（半正定値順序で）になることが証明されました。この関係は後の数値実験で実際に確認します。

バックワードパスのアルゴリズム

初期条件として \(\hat{x}_{T|T}\) と \(P_{T|T}\) （フォワードパスの最終結果）を設定し、\(k = T-1, T-2, \ldots, 0\) について式 \((8)\) 〜\((10)\) を逆順に適用します。

Python実装

1次元の位置・速度追跡問題を例として、カルマンフィルタとRTSスムーザを実装します。

状態空間モデルの定義

等速直線運動モデルを使用します。状態ベクトルは \(x = [位置, 速度]^T\) です。

import numpy as np
import matplotlib.pyplot as plt

# ---- 状態空間モデルの定義 ----
dt = 1.0  # サンプリング間隔

# 状態遷移行列（等速直線運動モデル）
F = np.array([[1, dt],
              [0, 1]])

# 観測行列（位置のみ観測）
H = np.array([[1, 0]])

# プロセスノイズ共分散
q = 0.1
Q = q * np.array([[dt**3/3, dt**2/2],
                   [dt**2/2, dt]])

# 観測ノイズ共分散
R = np.array([[1.0]])

n = 2  # 状態次元
m = 1  # 観測次元

カルマンフィルタ（フォワードパス）

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] = x0
    P_filt[0] = 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)
        innovation = y_obs[k] - H @ x_pred[k]
        x_filt[k + 1] = x_pred[k] + K @ innovation
        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] = x_filt[T]
    P_smooth[T] = P_filt[T]

    # バックワード再帰
    for k in range(T - 1, -1, -1):
        # スムーザゲイン（式8）
        G = P_filt[k] @ F.T @ np.linalg.inv(P_pred[k])

        # 平滑化状態推定（式9）
        x_smooth[k] = x_filt[k] + G @ (x_smooth[k + 1] - x_pred[k])

        # 平滑化共分散（式10）
        P_smooth[k] = P_filt[k] + G @ (P_smooth[k + 1] - P_pred[k]) @ G.T

    return x_smooth, P_smooth

シミュレーション

np.random.seed(42)
T = 50  # タイムステップ数

# 真の初期状態
x_true_init = np.array([0.0, 1.0])  # 位置0, 速度1

# 真の状態の生成
true_states = np.zeros((T + 1, n))
true_states[0] = x_true_init
for k in range(T):
    true_states[k + 1] = F @ true_states[k] + \
        np.random.multivariate_normal(np.zeros(n), Q)

# 観測の生成
y_obs = np.zeros((T, m))
for k in range(T):
    y_obs[k] = H @ true_states[k + 1] + \
        np.random.multivariate_normal(np.zeros(m), R)

# フィルタの初期推定
x0 = np.array([0.0, 0.0])
P0 = np.diag([1.0, 1.0])

# カルマンフィルタの実行
x_filt, P_filt, x_pred, P_pred = kalman_filter(y_obs, F, H, Q, R, x0, P0)

# RTSスムーザの実行
x_smooth, P_smooth = rts_smoother(x_filt, P_filt, x_pred, P_pred, F)

フィルタとスムーザの比較

time = np.arange(T + 1)

# ---- 位置の推定結果 ----
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(time, true_states[:, 0], "b-", linewidth=2, label="True state")
plt.scatter(np.arange(1, T + 1), y_obs[:, 0],
            c="gray", s=15, alpha=0.5, label="Measurements")
plt.plot(time, x_filt[:, 0], "r--", linewidth=1.5, label="KF estimate")
plt.plot(time, x_smooth[:, 0], "g-", linewidth=1.5, label="RTS estimate")
plt.xlabel("Time step")
plt.ylabel("Position")
plt.title("Position Estimation: Kalman Filter vs RTS Smoother")
plt.legend()
plt.grid(True)

# ---- 位置推定の誤差共分散 ----
plt.subplot(2, 1, 2)
kf_pos_var = np.array([P_filt[k, 0, 0] for k in range(T + 1)])
rts_pos_var = np.array([P_smooth[k, 0, 0] for k in range(T + 1)])
plt.plot(time, kf_pos_var, "r--", linewidth=1.5, label="KF variance")
plt.plot(time, rts_pos_var, "g-", linewidth=1.5, label="RTS variance")
plt.xlabel("Time step")
plt.ylabel("Position variance")
plt.title("Estimation Error Covariance Comparison")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.savefig("rts_smoother_result.png", dpi=150)
plt.show()

速度（直接観測されない状態）の推定結果を以下に示します。スムーザは観測されない状態成分に対しても大幅な精度改善を実現します。

RMSE比較

# 位置のRMSE
rmse_kf = np.sqrt(np.mean((true_states[1:, 0] - x_filt[1:, 0])**2))
rmse_rts = np.sqrt(np.mean((true_states[1:, 0] - x_smooth[1:, 0])**2))

# 速度のRMSE
rmse_kf_vel = np.sqrt(np.mean((true_states[1:, 1] - x_filt[1:, 1])**2))
rmse_rts_vel = np.sqrt(np.mean((true_states[1:, 1] - x_smooth[1:, 1])**2))

print(f"Position RMSE - KF: {rmse_kf:.4f}, RTS: {rmse_rts:.4f}")
print(f"Velocity RMSE - KF: {rmse_kf_vel:.4f}, RTS: {rmse_rts_vel:.4f}")
print(f"Position improvement: {(1 - rmse_rts / rmse_kf) * 100:.1f}%")
print(f"Velocity improvement: {(1 - rmse_rts_vel / rmse_kf_vel) * 100:.1f}%")

RTSスムーザの推定は、カルマンフィルタの推定よりも真の状態に近く、特にデータの端部以外では大幅な精度改善が確認できます。誤差共分散のプロットからも、スムーザの不確実性が常にフィルタ以下であることが分かります。

数値検証：導出した性質を実際に確かめる

前節で導出した「\(P_{k|T} \leq P_{k|k}\) が常に成り立つ」という性質と、区間端点での挙動を、上のシミュレーションデータを使って実際に検証します。

\(P_{k|T} \leq P_{k|k}\) の検証（半正定値順序）

# 全時刻で P_filt[k] - P_smooth[k] が半正定値（固有値がすべて非負）であることを確認
eigmins = np.array([
    np.linalg.eigvalsh(P_filt[k] - P_smooth[k]).min() for k in range(T + 1)
])
print(f"min eigenvalue of P_filt - P_smooth over all k: {eigmins.min():.6e}")
print(f"any negative below -1e-10 (floating point tolerance)?: {np.any(eigmins < -1e-10)}")

# 区間端点・中央での分散比（スムーザ/フィルタ）
for k in [0, T // 2, T]:
    ratio = P_smooth[k, 0, 0] / P_filt[k, 0, 0]
    print(f"k={k:2d}: position variance ratio (smooth/filt) = {ratio:.4f}")

実行結果：

min eigenvalue of P_filt - P_smooth over all k: -2.701731e-17
any negative below -1e-10 (floating point tolerance)?: False
k= 0: position variance ratio (smooth/filt) = 0.5112
k=25: position variance ratio (smooth/filt) = 0.3624
k=50: position variance ratio (smooth/filt) = 1.0000

最小固有値は \(-2.7 \times 10^{-17}\) と、倍精度浮動小数点の丸め誤差の範囲内でゼロに一致しており、全時刻で半正定値性（証明した \(P_{k|T} \leq P_{k|k}\) ）が数値的にも確認できました。分散比に注目すると、区間の端点 \(k=T\) ではちょうど \(1.0000\) （等号成立）になっています。これは証明の base case（\(P_{T|T}=P_{T|T}\) ）が示す通り、最終時刻には未来の観測が存在しないためスムーザがフィルタと完全に一致するという、区間端点特有のエッジケースです。一方、区間内部（\(k=25\) ）では分散が \(36\%\) まで縮小し、始点付近（\(k=0\) ）でも約半分に縮小しています。始点でスムーザの改善効果が中間点よりやや小さいのは、初期共分散 \(P_0\) （本記事では \(\mathrm{diag}(1,1)\) ）の不確実性が残るためです。

エッジケース：観測欠損時のスムーザ

実運用では、通信断やセンサ故障により、一部の時刻で観測が得られないことがあります。カルマンフィルタは観測がない時刻について単純に更新ステップをスキップし（\(\hat{x}_{k|k} = \hat{x}_{k|k-1}\) 、\(P_{k|k} = P_{k|k-1}\) ）予測のみで代替しますが、この間はフィルタの不確実性が単調に増大し続けます。RTSスムーザは後ろ向き再帰の式 \((8)\) 〜\((10)\) をそのまま適用できるため（観測の有無は \(x_{k|k}, P_{k|k}, x_{k|k-1}, P_{k|k-1}\) の値に反映済みで、再帰式自体は変更不要）、欠損区間の前後の情報を使って大幅に精度を回復できます。

def kalman_filter_missing(y_obs, mask, F, H, Q, R, x0, P0):
    """観測欠損対応カルマンフィルタ。mask[k]=False の時刻は更新をスキップする"""
    T = len(y_obs)
    n = len(x0)
    x_filt = np.zeros((T + 1, n))
    P_filt = np.zeros((T + 1, n, n))
    x_pred = np.zeros((T, n))
    P_pred = np.zeros((T, n, n))
    x_filt[0], P_filt[0] = x0, P0
    for k in range(T):
        x_pred[k] = F @ x_filt[k]
        P_pred[k] = F @ P_filt[k] @ F.T + Q
        if mask[k]:
            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]
        else:
            # 観測なし: 予測をそのままフィルタ推定として採用（不確実性は増大したまま）
            x_filt[k + 1] = x_pred[k]
            P_filt[k + 1] = P_pred[k]
    return x_filt, P_filt, x_pred, P_pred

# 時刻20〜29（10ステップ連続、全50ステップ中20%）の観測が欠損したとする
gap_start, gap_end = 20, 30
mask = np.ones(T, dtype=bool)
mask[gap_start:gap_end] = False
gap_idx = np.arange(gap_start + 1, gap_end + 1)  # 対応する状態インデックス

x_filt_m, P_filt_m, x_pred_m, P_pred_m = kalman_filter_missing(y_obs, mask, F, H, Q, R, x0, P0)
x_smooth_m, P_smooth_m = rts_smoother(x_filt_m, P_filt_m, x_pred_m, P_pred_m, F)

rmse_filt_gap = np.sqrt(np.mean((true_states[gap_idx, 0] - x_filt_m[gap_idx, 0]) ** 2))
rmse_smooth_gap = np.sqrt(np.mean((true_states[gap_idx, 0] - x_smooth_m[gap_idx, 0]) ** 2))
print(f"Position RMSE in gap - filter: {rmse_filt_gap:.4f}, smoother: {rmse_smooth_gap:.4f}")
print(f"Position variance at gap middle (k=25) - filter: {P_filt_m[25,0,0]:.4f}, smoother: {P_smooth_m[25,0,0]:.4f}")
print(f"Position variance at gap end (k=30)    - filter: {P_filt_m[30,0,0]:.4f}, smoother: {P_smooth_m[30,0,0]:.4f}")

実行結果（この1系列では）：

Position RMSE in gap - filter: 0.6071, smoother: 0.5778
Position variance at gap middle (k=25) - filter: 12.0439, smoother: 1.9103
Position variance at gap end (k=30)    - filter: 58.9471, smoother: 0.6918

分散に注目すると、欠損区間の終端（\(k=30\) ）でフィルタの位置分散が \(58.9\) まで膨れ上がっているのに対し、スムーザは \(0.69\) と平常時とほぼ同水準に保たれています。これは、スムーザが欠損区間の後にある観測（\(k=30\) 以降）から逆向きに情報を伝播できるためです。ただし、この1系列だけではRMSEの改善が\(4.8\%\) と小さく、たまたま速度成分ではフィルタの方がわずかに良い（サンプリングの偶然による揺らぎ）という結果も出ました。これは欠損区間内のサンプル数が10点と少ないため、単一系列のRMSEは統計的に不安定であることを示しています。そこで、200試行のモンテカルロ平均で傾向を確認します。

rng = np.random.default_rng(123)
n_trials = 200
pos_f, pos_s, vel_f, vel_s = [], [], [], []
for trial in range(n_trials):
    ts = np.zeros((T + 1, n))
    ts[0] = [0.0, 1.0]
    for k in range(T):
        ts[k + 1] = F @ ts[k] + rng.multivariate_normal(np.zeros(n), Q)
    yo = np.zeros((T, m))
    for k in range(T):
        yo[k] = H @ ts[k + 1] + rng.multivariate_normal(np.zeros(m), R)

    xf, Pf, xp, Pp = kalman_filter_missing(yo, mask, F, H, Q, R, x0, P0)
    xs, Ps = rts_smoother(xf, Pf, xp, Pp, F)

    pos_f.append(np.sqrt(np.mean((ts[gap_idx, 0] - xf[gap_idx, 0]) ** 2)))
    pos_s.append(np.sqrt(np.mean((ts[gap_idx, 0] - xs[gap_idx, 0]) ** 2)))
    vel_f.append(np.sqrt(np.mean((ts[gap_idx, 1] - xf[gap_idx, 1]) ** 2)))
    vel_s.append(np.sqrt(np.mean((ts[gap_idx, 1] - xs[gap_idx, 1]) ** 2)))

pos_f, pos_s, vel_f, vel_s = map(np.array, (pos_f, pos_s, vel_f, vel_s))
print(f"Position RMSE in gap - filter mean: {pos_f.mean():.4f}, smoother mean: {pos_s.mean():.4f}")
print(f"Velocity RMSE in gap - filter mean: {vel_f.mean():.4f}, smoother mean: {vel_s.mean():.4f}")
print(f"fraction of trials smoother beats filter (position): {(pos_s < pos_f).mean()*100:.1f}%")
print(f"fraction of trials smoother beats filter (velocity): {(vel_s < vel_f).mean()*100:.1f}%")

実行結果：

Position RMSE in gap - filter mean: 3.6140, smoother mean: 1.0532
Velocity RMSE in gap - filter mean: 0.7126, smoother mean: 0.3244
fraction of trials smoother beats filter (position): 87.5%
fraction of trials smoother beats filter (velocity): 92.0%

200試行平均では、欠損区間の位置RMSEが \(70.9\%\) 、速度RMSEが \(54.5\%\) 改善し、スムーザがフィルタに勝つ試行の割合も位置で\(87.5\%\) 、速度で\(92.0\%\) に達しました。単一系列では偶然フィルタが上回ることもありますが（速度成分で実際に発生しました）、平均を取れば理論通りスムーザが一貫して優位であることが確認できます。単一の実行結果だけで手法の優劣を判断せず、複数試行の平均で評価すべきという一般的な教訓が、ここでも当てはまります。

下図は、この欠損シナリオでの位置推定と \(\pm 2\sigma\) 不確実性帯を示したものです。欠損区間（灰色）でフィルタの不確実性帯（赤）が大きく広がる一方、スムーザの不確実性帯（緑）はほぼ変化しないことが視覚的に確認できます。

数値安定性：悪条件な場合とJoseph形式

バックワード再帰の式 \((8)\) は \(P_{k+1|k}^{-1}\) の計算を要求します。\(P_{k+1|k}\) が悪条件（ほぼ特異）になると、この逆行列計算で数値誤差が大きく増幅されるおそれがあります。

いつ悪条件になるか

プロセスノイズ \(Q\) が観測ノイズ \(R\) に対して極端に小さい場合（ほぼ決定論的なシステム、あるいはある状態成分がノイズをほとんど受けない場合）、\(P_{k+1|k} = F P_{k|k} F^T + Q\) の固有値のスケールが大きく乖離し、条件数 \(\mathrm{cond}(P_{k+1|k}) = \lambda_{\max}/\lambda_{\min}\) が急増します。逆行列・連立方程式の相対誤差は一般に機械イプシロン \(\epsilon\) と条件数の積 \(O(\epsilon \cdot \mathrm{cond}(P))\) のオーダーでスケールするため（数値線形代数の標準的な誤差解析、Golub & Van Loan (2013) 等）、条件数が極端に大きいケースでは有効桁が大きく失われます。次の実験で、条件数と実際の相対誤差の関係を人工的な悪条件行列で直接確認します。

rng2 = np.random.default_rng(0)

def make_ill_conditioned(cond_number, size=2):
    """条件数を指定した対称正定値行列を人工的に構成する（直交行列由来なので厳密逆行列が既知）"""
    A = rng2.standard_normal((size, size))
    Qm, _ = np.linalg.qr(A)
    eigs = np.geomspace(1.0, cond_number, size)
    P = Qm @ np.diag(eigs) @ Qm.T
    P_inv_exact = Qm @ np.diag(1.0 / eigs) @ Qm.T  # 直交行列なので厳密逆行列
    return P, P_inv_exact

print(f"{'cond(P)':>10s} | {'rel. error inv()':>18s} | {'rel. error solve()':>18s}")
for cond_target in [1e2, 1e4, 1e6, 1e8, 1e10, 1e12, 1e14]:
    P, P_inv_exact = make_ill_conditioned(cond_target)
    B = rng2.standard_normal((2, 2))       # P_{k|k} F^T に相当するプレースホルダ
    G_ref = B @ P_inv_exact                # 厳密解

    G_inv = B @ np.linalg.inv(P)           # 素朴な明示的逆行列（式8の実装通り）
    G_solve = np.linalg.solve(P.T, B.T).T  # 連立方程式として解く（逆行列を作らない）

    err_inv = np.linalg.norm(G_inv - G_ref) / np.linalg.norm(G_ref)
    err_solve = np.linalg.norm(G_solve - G_ref) / np.linalg.norm(G_ref)
    print(f"{np.linalg.cond(P):10.3e} | {err_inv:18.3e} | {err_solve:18.3e}")

実行結果：

   cond(P) |   rel. error inv() | rel. error solve()
 1.000e+02 |          8.239e-16 |          6.652e-16
 1.000e+04 |          1.181e-13 |          1.181e-13
 1.000e+06 |          2.496e-11 |          3.042e-12
 1.000e+08 |          1.404e-09 |          1.404e-09
 1.000e+10 |          5.397e-07 |          2.321e-07
 1.000e+12 |          1.288e-05 |          1.288e-05
 1.001e+14 |          2.496e-04 |          2.496e-04

条件数が4桁増えるごとに相対誤差もおおむね4桁増えており、\(O(\epsilon \cdot \mathrm{cond}(P))\) という理論的なスケーリングとほぼ整合しています。np.linalg.solve は明示的な逆行列を作らずに連立方程式として解くため、条件数 \(10^6\) 〜\(10^{10}\) の領域ではやや誤差が小さくなる場合がありますが（例えば \(10^6\) で約8倍の改善）、根本的な悪条件性そのものを解消するわけではありません。

Joseph形式による対称性の保護

フィルタの更新式 \((7)\) の \(P_{k|k} = (I-K_kH)P_{k|k-1}\) は減算を含むため、丸め誤差によって対称性が壊れることがあります。これに対しJoseph形式

は2つの合同変換の和として書かれているため、\(P_{k|k-1}\) と \(R\) が半正定値である限り、丸め誤差があっても構造的に対称性・半正定値性が保たれやすいという利点があります（計算コストは通常形の約2倍）。実際に、単精度浮動小数点（float32）・\(T=2000\) ステップ・プロセスノイズ約\(10^{-14}\) （ほぼ決定論的な等速直線運動）・観測ノイズ\(R=10^{-6}\) （非常に精密な観測）という悪条件な設定でフィルタを実行し、共分散の非対称性 \(\|P-P^\top\|\) を比較すると、通常形の最大非対称性が \(2.090\times10^{-12}\) であったのに対し、Joseph形式では \(1.608\times10^{-13}\) と約13倍小さくなりました。一方で、この間の予測共分散 \(P_{k+1|k}\) の条件数は両形式とも最大で約 \(2\times10^6\) に達しており、Joseph形式は対称性の保護には有効ですが、\(P_{k+1|k}^{-1}\) の悪条件性そのものは解消しないことも確認されました。

実務上の指針としては、(1) プロセスノイズに下限を設けて（\(Q \leftarrow Q + \delta I\) 、\(\delta\) は微小値）\(P_{k+1|k}\) が特異に近づくのを防ぐ、(2) 明示的な np.linalg.inv の代わりに np.linalg.solve（またはCholesky分解ベースの解法）を使う、(3) 更新のたびに \(P \leftarrow (P+P^\top)/2\) で明示的に対称化する、(4) フィルタの共分散更新にはJoseph形式を用いる、という4点を組み合わせるのが安全です。

近年の研究動向

RTSスムーザは1965年の原論文以来ほぼ形を変えずに使われ続けていますが、2020年代に入り2つの方向で発展が続いています。

非線形系への拡張：本記事の\(G_k = P_{k|k}F^\top P_{k+1|k}^{-1}\) という漸化式は線形ガウスモデルを前提としていますが、遷移関数・観測関数が非線形な場合には Extended RTS Smoother（EKS、ヤコビアンで線形近似）や Unscented RTS Smoother（URTS、シグマポイントで非線形変換の統計量を近似、Särkkä, 2008）への拡張が知られています。これらの導出・実装・強非線形ベンチマークでの精度比較は、 非線形カルマンスムーザ（EKS/URTS）の理論とPython実装 で詳しく扱っているため、本記事では概要にとどめます。

深層学習との融合：状態遷移・観測モデルが部分的にしか分からない（model mismatch）実問題に対し、古典的なRTS漸化式の構造を保ったまま、スムーザゲインの計算をニューラルネットワークで補完する RTSNet（Revach, Ni, Shlezinger, van Sloun, & Eldar, 2023, deep unfolding によりRTS再帰へ学習可能なモジュールを組み込む手法。姉妹モデルの KalmanNet, Revach et al., 2022 のスムーザ版）が提案されています。RTSNetは、モデルが一部不明・非線形な設定でも、古典的スムーザの計算効率と解釈可能性を保持しながら、モデル誤差を学習によって補正できる点が特徴です。同様の設計思想は Latent-KalmanNet（高次元信号からの学習型フィルタリング、Revach et al., 2023）などにも受け継がれており、「モデルベースの再帰構造にデータ駆動の補正項を組み込む」というハイブリッドアプローチが、深層状態空間モデル（Deep State Space Model）研究の一潮流として定着しつつあります。本記事で導出した後ろ向き再帰の数式的な骨格は、こうした学習型手法においても基盤として使われ続けています。

まとめ

RTSスムーザは、カルマンフィルタのフォワードパスに対してバックワードパスを追加するだけで実現でき、実装コストが低い割に推定精度を大きく向上させます。リアルタイム性が不要なバッチ処理のシナリオでは、まずRTSスムーザの適用を検討すべきです。

• 後ろ向き再帰が正しい平滑化分布を与えることを、ベイズの定理とマルコフ性（\(p(x_k\mid x_{k+1},z_{1:T})=p(x_k\mid x_{k+1},z_{1:k})\) ）から厳密に証明し、平滑化ゲイン \(G_k\) をガウス分布の条件付け公式から導出した
• \(P_{k|T}\leq P_{k|k}\) を後ろ向き帰納法で証明し、数値実験でも最小固有値が\(-2.7\times10^{-17}\) （浮動小数点誤差の範囲内でゼロ）であることを確認した
• 観測欠損（10ステップの欠損区間）でも、200試行モンテカルロ平均でスムーザが位置RMSEを\(70.9\%\) 、速度RMSEを\(54.5\%\) 改善することを確認した
• プロセスノイズが極端に小さい悪条件なケースでは、\(P_{k+1|k}\) の条件数が\(10^6\) 規模に達し、Joseph形式が共分散の非対称性を約13倍抑制する一方、悪条件性そのものは解消しないことを確認した

なお、本記事では線形モデルに限定しましたが、非線形モデルに対しては Extended Kalman Smoother や Unscented Kalman Smoother、 CKF ベースのスムーザなどの拡張も存在します（詳細は上記「近年の研究動向」を参照）。

おすすめ書籍

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

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

関連記事

• カルマンスムーザの比較 - RTSスムーザを含む複数のスムージング手法を横断比較しています。本記事の次に読むと手法選択の指針が得られます。
• 非線形カルマンスムーザ（EKS/URTS）の理論とPython実装 - 本記事の後ろ向き再帰を非線形モデルへ拡張したExtended RTS SmootherとUnscented RTS Smootherを導出・実装し、強非線形ベンチマークで精度を比較しています。
• カルマンフィルタの理論とPython実装 - RTSスムーザの基盤となるカルマンフィルタの理論と実装を解説しています。
• 拡張カルマンフィルタ（EKF）の理論とPython実装 - 非線形モデルに対応するEKFを解説しています。非線形スムーザの基盤にもなります。
• Unscented Kalman Filter（UKF）の理論とPython実装 - UKFベースのスムーザへの拡張の基盤となるUKFを解説しています。
• 粒子フィルタのPython実装：リサンプリング手法の比較 - 粒子スムーザの基盤となる粒子フィルタを解説しています。
• 時系列データの異常検知：統計的手法からカルマンフィルタまで - RTSスムーザはオフライン異常検知の前処理としても有用です。
• ウィーナーフィルタの理論とPython実装 - 周波数領域での定常過程の最適フィルタ。RTSの時間領域・状態空間ベースのスムーザと数学的対応関係があります。
• ガウス過程回帰の理論とPython実装 - 状態空間モデルとガウス過程の双対性が知られており、RTSスムーザはGP回帰のオンライン版と見ることができます。
• サビツキー・ゴーレイフィルタの理論とPython実装 - オフラインでの平滑化手法。状態空間モデルを使わない局所多項式フィッティングとの対比が学べます。

参考文献

• Rauch, H. E., Tung, F., & Striebel, C. T. (1965). “Maximum likelihood estimates of linear dynamic systems.” AIAA Journal, 3(8), 1445-1450.
• Särkkä, S. (2013). Bayesian Filtering and Smoothing. Cambridge University Press.
• Särkkä, S. (2008). “Unscented Rauch-Tung-Striebel Smoother.” IEEE Transactions on Automatic Control, 53(3), 845-849.
• Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
• Revach, G., Ni, X., Shlezinger, N., van Sloun, R. J. G., & Eldar, Y. C. (2023). “RTSNet: Learning to Smooth in Partially Known State-Space Models.” IEEE Transactions on Signal Processing, 71, 4441-4456.