RLS(逐次最小二乗法)アルゴリズムの理論とPython実装:カルマンフィルタとの等価性と忘却係数

RLS適応フィルタを行列反転補題(Sherman-Morrison公式)からnumpy.outer・numpy.eyeでフルスクラッチ実装。200試行モンテカルロでLMSとの収束速度を比較し、忘却係数λ=1のRLSとカルマンフィルタが数値的に完全一致することを実測、非定常環境での忘却係数λの追従性と定常誤差のトレードオフも検証します。

はじめに

https://yuhi-sa.github.io/posts/20260310_adaptive_filter/1/ と https://yuhi-sa.github.io/posts/20260501_lms_nlms/1/ では、LMS・NLMSアルゴリズムを中心に適応フィルタを解説しました。両記事とも最後に「RLSはLMSより高速に収束するが計算コストが高い」「RLSはカルマンフィルタと数学的に等価」と触れるにとどめ、導出には踏み込んでいませんでした。本記事ではその続きとして、**RLS(Recursive Least Squares、逐次最小二乗法)**を行列反転補題から導出し、Pythonでフルスクラッチ実装します。さらに、RLSとカルマンフィルタが特定の条件下で数値的に完全一致することを実験で確かめ、忘却係数 \(\lambda\) が定常誤差と追従性のトレードオフをどう支配するかを検証します。

問題設定:指数重み付き最小二乗

RLSは、LMSのような瞬時勾配の確率的近似ではなく、過去の全観測に対する重み付き二乗誤差を厳密に最小化します。長さ \(M\) のFIRフィルタ係数 \(\mathbf{w}\) に対し、時刻 \(n\) でのコスト関数を

\[ J(\mathbf{w}, n) = \sum_{i=1}^{n} \lambda^{n-i} \bigl[d(i) - \mathbf{x}(i)^\top \mathbf{w}\bigr]^2 \tag{1} \]

と定義します。\(\lambda \in (0, 1]\) は忘却係数(forgetting factor)で、過去のサンプルの寄与を指数的に減衰させます。\(\lambda = 1\) ならすべてのサンプルを等しく扱う通常の最小二乗、\(\lambda < 1\) なら直近のサンプルを重視する適応的な推定になります。

\(J(\mathbf{w}, n)\) を \(\mathbf{w}\) で微分してゼロと置くと、正規方程式

\[ \mathbf{R}(n) \mathbf{w}(n) = \mathbf{r}(n), \qquad \mathbf{R}(n) = \sum_{i=1}^n \lambda^{n-i} \mathbf{x}(i)\mathbf{x}(i)^\top, \quad \mathbf{r}(n) = \sum_{i=1}^n \lambda^{n-i} \mathbf{x}(i) d(i) \tag{2} \]

が得られます。\(\mathbf{R}(n)\) を毎回逆行列計算すると \(O(M^3)\) かかりますが、次の行列反転補題を使えば \(O(M^2)\) の再帰更新に落とし込めます。

行列反転補題によるO(M²)再帰更新

\(\mathbf{R}(n)\) は前時刻から次のように更新できます。

\[ \mathbf{R}(n) = \lambda \mathbf{R}(n-1) + \mathbf{x}(n)\mathbf{x}(n)^\top \tag{3} \]

\(P(n) = \mathbf{R}(n)^{-1}\) とおき、Sherman-Morrison公式(行列反転補題)

\[ (A + \mathbf{u}\mathbf{v}^\top)^{-1} = A^{-1} - \frac{A^{-1}\mathbf{u}\mathbf{v}^\top A^{-1}}{1 + \mathbf{v}^\top A^{-1}\mathbf{u}} \tag{4} \]

式(4)の検証。 この式は天下り的に見えますが、右辺が本当に左辺の逆行列になっていることは直接代入で確認できます。\(s = 1 + \mathbf{v}^\top A^{-1}\mathbf{u}\) (スカラー)とおき、式(4)の右辺を \(B\) と書くと、

\[ (A + \mathbf{u}\mathbf{v}^\top) B = A A^{-1} - \frac{A A^{-1}\mathbf{u}\mathbf{v}^\top A^{-1}}{s} + \mathbf{u}\mathbf{v}^\top A^{-1} - \frac{\mathbf{u}\mathbf{v}^\top A^{-1}\mathbf{u}\mathbf{v}^\top A^{-1}}{s} \]

\(\mathbf{v}^\top A^{-1}\mathbf{u}\) はスカラーなのでベクトルの外に出せて、\(\mathbf{v}^\top A^{-1}\mathbf{u} = s - 1\) を使うと、

\[ (A + \mathbf{u}\mathbf{v}^\top) B = I - \frac{\mathbf{u}\mathbf{v}^\top A^{-1}}{s} + \mathbf{u}\mathbf{v}^\top A^{-1} - \frac{(s-1)\,\mathbf{u}\mathbf{v}^\top A^{-1}}{s} = I + \mathbf{u}\mathbf{v}^\top A^{-1}\cdot\frac{-1 + s - (s-1)}{s} = I \]

最後の括弧内の分子がちょうど \(0\) になるため、\((A+\mathbf{u}\mathbf{v}^\top)B = I\) が厳密に成り立ちます(\(B(A+\mathbf{u}\mathbf{v}^\top) = I\) も同様の計算で示せます)。この直接検証により、以下でSherman-Morrison公式を適用して得られる再帰更新式(5)〜(7)が近似ではなく厳密な等式であることが保証されます。

式(4)を \(A = \lambda \mathbf{R}(n-1)\) 、\(\mathbf{u} = \mathbf{v} = \mathbf{x}(n)\) に適用すると、

\[ P(n) = \frac{1}{\lambda}\left[P(n-1) - \frac{P(n-1)\mathbf{x}(n)\mathbf{x}(n)^\top P(n-1)}{\lambda + \mathbf{x}(n)^\top P(n-1)\mathbf{x}(n)}\right] \tag{5} \]

が得られます。逆行列を直接計算せず、既知の \(P(n-1)\) から行列積とスカラー除算だけで \(P(n)\) を求められるのがポイントです。ここでカルマンゲインにあたる量

\[ \mathbf{k}(n) = \frac{P(n-1)\mathbf{x}(n)}{\lambda + \mathbf{x}(n)^\top P(n-1)\mathbf{x}(n)} \tag{6} \]

を定義すると、式(5)は \(P(n) = \frac{1}{\lambda}\bigl[P(n-1) - \mathbf{k}(n)\mathbf{x}(n)^\top P(n-1)\bigr]\) と簡潔に書けます。重み更新は

\[ e(n) = d(n) - \mathbf{x}(n)^\top \mathbf{w}(n-1), \qquad \mathbf{w}(n) = \mathbf{w}(n-1) + \mathbf{k}(n) e(n) \tag{7} \]

です。1ステップあたりの演算量は \(\mathbf{k}(n)\) の計算が支配的で \(O(M^2)\) 、LMS/NLMSの \(O(M)\) より重くなります。

Python実装

import numpy as np

def rls_filter(xs, d, M, lam=1.0, delta=1.0):
    """RLSアルゴリズムによる適応フィルタ

    Parameters
    ----------
    xs : ndarray, shape (N, M)
        各時刻の入力ベクトル
    d : ndarray, shape (N,)
        所望信号
    M : int
        フィルタ長
    lam : float
        忘却係数(0 < lam <= 1)
    delta : float
        初期共分散の正則化パラメータ(P(0) = I / delta)

    Returns
    -------
    W : ndarray, shape (N, M)
        各時刻のフィルタ係数の推移
    """
    w = np.zeros(M)
    P = np.eye(M) / delta
    W = np.zeros((len(d), M))
    for n in range(len(d)):
        x = xs[n]
        Px = P @ x
        k = Px / (lam + x @ Px)
        e = d[n] - x @ w
        w = w + k * e
        P = (P - np.outer(k, x @ P)) / lam
        W[n] = w
    return W

np.outer(k, x @ P) が式(5)のランク1更新項、Px / (lam + x @ Px) が式(6)のゲインに対応します。

数値実験1:LMSとの収束速度比較

\(M=8\) タップのFIRシステム同定問題(SNR 20dB)で、LMS(ステップサイズ \(\mu=0.05\) )とRLS(\(\lambda=1\) )を200試行モンテカルロ平均した学習曲線を比較しました。乱数シードを固定した完全な実験コードは以下の通りです(コピーしてそのまま実行すれば同じ数値が再現できます)。

np.random.seed(42)
M, N, n_trials, snr_db = 8, 150, 200, 20

true_w = np.random.randn(M)
true_w /= np.linalg.norm(true_w)
noise_power = 1.0 / (10 ** (snr_db / 10))  # 信号パワーを1に正規化

lms_sq_err = np.zeros(N)
rls_sq_err = np.zeros(N)
for trial in range(n_trials):
    xs = np.random.randn(N, M)
    d = xs @ true_w + np.random.randn(N) * np.sqrt(noise_power)
    W_lms = lms_filter(xs, d, M, mu=0.05)
    W_rls = rls_filter(xs, d, M, lam=1.0, delta=1.0)
    lms_sq_err += np.sum((W_lms - true_w) ** 2, axis=1)
    rls_sq_err += np.sum((W_rls - true_w) ** 2, axis=1)

lms_mse = lms_sq_err / n_trials
rls_mse = rls_sq_err / n_trials

lms_filter はLMSの標準実装 w += mu * e * x です。)実行結果は次の通りです。

サンプル数 \(n\)LMSのMSERLSのMSE
200.20530.0148
500.02410.0025
1000.00310.0010

雑音床の2倍(MSE ≈ 0.02)に到達するまでのサンプル数は、LMSが54サンプル、RLSが18サンプルでした。RLSはLMSの3分の1程度のサンプル数で収束しており、理論通り「LMSは線形収束、RLSは二次収束(入力の統計に依存せず最短時間で正規方程式の解に到達)」という特性が数値的に確認できます。一方でRLSは1ステップ \(O(M^2)\) のためLMSの \(O(M)\) より計算コストが高く、この収束速度と計算量のトレードオフが用途選択の基準になります。下図左は、この学習曲線をそのまま可視化したものです(右パネルは次節の忘却係数の実験に対応します)。

RLSとLMSの収束速度比較、および忘却係数λの効果

RLSとカルマンフィルタの等価性

https://yuhi-sa.github.io/posts/20260224_kalman_filter/1/ で導出したカルマンフィルタの更新式(予測・イノベーション・ゲイン・更新)を思い出してください。RLSは、以下のようにカルマンフィルタの特殊ケースとして定式化できます。

状態空間モデルとして見たRLS:

カルマンフィルタの記号RLSでの対応意味
状態 \(\mathbf{x}_k\)フィルタ係数 \(\mathbf{w}\)推定したいパラメータ
状態遷移 \(A\)\(I\) (単位行列)係数は時間的に一定と仮定
プロセスノイズ \(Q\)\(0\)係数の真値は変化しない
観測行列 \(H\)\(\mathbf{x}(n)^\top\)時刻ごとに変わる入力ベクトル
観測ノイズ \(R\)\(\lambda\) (定数)忘却係数がそのまま観測分散に対応
初期共分散 \(P_0\)\(I / \delta\)RLSの正則化パラメータと一致

このマッピングでカルマンフィルタの更新式(\(K_k = P_{k|k-1}H^\top S_k^{-1}\) 等)を書き下すと、式(6)(7)のRLS更新式と完全に一致します。実際に両者を実装して同じデータに適用し、係数の推移を比較しました。

def kalman_param_filter(xs, d, M, R=1.0, delta=1.0):
    w = np.zeros(M)
    P = np.eye(M) / delta
    W = np.zeros((len(d), M))
    for n in range(len(d)):
        x = xs[n]
        Px = P @ x
        S = x @ Px + R          # イノベーション共分散
        K = Px / S               # カルマンゲイン
        e = d[n] - x @ w
        w = w + K * e
        P = P - np.outer(K, x @ P)
        W[n] = w
    return W

\(M=4\) タップ、500サンプルのシステム同定問題で \(\lambda = 1\) 、\(R = 1\) 、\(\delta\) を揃え、np.random.seed(0) で乱数を固定して両者を実行しました。コピーしてそのまま実行すれば同じ数値が再現できる完全なコードは次の通りです。

np.random.seed(0)
M, N = 4, 500
true_w = np.array([1.0, -0.5, 0.3, 0.2])
xs = np.random.randn(N, M)
d = xs @ true_w + np.random.randn(N) * 1.0  # R=1 なので観測ノイズの標準偏差は1

W_rls = rls_filter(xs, d, M, lam=1.0, delta=1.0)
W_kf = kalman_param_filter(xs, d, M, R=1.0, delta=1.0)
print(f"max |RLS - KF| over all n, all taps (lambda=1, R=1): {np.max(np.abs(W_rls - W_kf))}")
print(f"final RLS weights: {W_rls[-1]}")
print(f"final KF  weights: {W_kf[-1]}")
print(f"true weights     : {true_w}")

実行結果:

max |RLS - KF| over all n, all taps (lambda=1, R=1): 0.0
final RLS weights: [ 0.9475194  -0.44072455  0.3474052   0.23557128]
final KF  weights: [ 0.9475194  -0.44072455  0.3474052   0.23557128]
true weights     : [ 1.  -0.5  0.3  0.2]

**全時刻・全タップで誤差が完全にゼロ(機械精度で一致)**となり、RLSがカルマンフィルタの特殊ケースであることが数値的に裏付けられました。この等価性は実務上も有用で、RLSの忘却係数を調整する代わりに、より一般的なカルマンフィルタの枠組み(プロセスノイズ \(Q\) を明示的に導入する)に切り替えれば、非定常性のモデル化を柔軟に行えます。

この等価性は2020年代に入っても発展を続けている研究テーマです。Lai and Bernstein (2024) は、様々な忘却拡張を持つRLS(generalized forgetting RLS)を統一的に再帰的最小化として捉える「Kalman filter least squares(KFLS)」コスト関数を提案し、本記事で数値的に確認した \(\lambda=1\) の場合だけでなく、より一般的な忘却拡張のRLS変種もすべてカルマンフィルタの特殊ケースとして導出できることを示しました。さらに、この統一的視点を利用して、間欠的で未モデル化の外乱(intermittent, unmodeled collisions)を受ける質量-バネ-ダンパ系に対し、従来のRLSより高精度な状態推定を達成する適応カルマンフィルタの新しいクラスを構築しています。本記事のRLS-カルマン等価性の数値実験は、この統一理論の最も単純な特殊ケース(\(Q=0\) 、定数の観測ノイズ)に相当します。

忘却係数λの効果:定常性 vs 追従性

\(\lambda < 1\) の近似的な直感は「プロセスノイズ \(Q > 0\) を許容するカルマンフィルタに近づく」というものです。これを検証するため、系の途中(1000サンプル目、\(N=1200\) )でフィルタ係数の真値が \([1.0, -0.5, 0.3, 0.2]\) から \([-0.8, 0.6, -0.2, 0.5]\) へ急変する非定常システムを用意し、\(\lambda\) を変えてRLSを適用しました。単一試行では観測ノイズの影響で誤差が非単調になりうるため、np.random.seed(1) を起点に200試行モンテカルロ平均を取っています。コピーしてそのまま実行すれば同じ数値が再現できる完全なコードは次の通りです。

np.random.seed(1)
M, N, jump_n, n_trials = 4, 1200, 1000, 200
true_w1 = np.array([1.0, -0.5, 0.3, 0.2])
true_w2 = np.array([-0.8, 0.6, -0.2, 0.5])
true_w_t = np.tile(true_w1, (N, 1))
true_w_t[jump_n:] = true_w2  # 1000サンプル目で真値が急変

lams = [1.000, 0.995, 0.980, 0.900]
steady_acc = {lam: 0.0 for lam in lams}
after_acc = {lam: 0.0 for lam in lams}
for trial in range(n_trials):
    xs = np.random.randn(N, M)
    d = np.sum(xs * true_w_t, axis=1) + np.random.randn(N) * 1.0
    for lam in lams:
        W = rls_filter(xs, d, M, lam=lam, delta=1.0)
        err = np.sum((W - true_w_t) ** 2, axis=1)  # ||w(n) - w*(n)||^2
        steady_acc[lam] += np.mean(err[jump_n - 50:jump_n])  # ジャンプ直前50サンプルの平均誤差
        after_acc[lam] += err[jump_n + 9]                     # ジャンプ後10サンプルの誤差

for lam in lams:
    print(f"lambda={lam}: steady={steady_acc[lam]/n_trials:.4f}  after10={after_acc[lam]/n_trials:.4f}")
\(\lambda\)ジャンプ直前の定常誤差ジャンプ後10サンプルの誤差
1.0000.00444.7116
0.9950.00984.4069
0.9800.03983.4944
0.9000.22201.3792

\(\lambda = 1\) は定常状態での誤差が最小(過去の全データを均等に使うため統計的に最も効率的)ですが、システムが変化した後の追従が最も遅く、ジャンプからわずか10サンプルの時点では新しい値へほとんど収束していません。逆に \(\lambda = 0.9\) は定常誤差が約50倍大きい(直近10サンプル程度しか実質的に使っていないのと同等)反面、変化への追従は最も速く、同じ10サンプル後の誤差はλ=1の場合の3分の1以下まで下がっています。\(\lambda\) は「定常時の精度」と「非定常環境への追従性」を直接トレードオフする設計パラメータであることが、この実験から定量的に確認できます。上図右パネルは、この4つの \(\lambda\) について係数誤差 \(\lVert \mathbf{w}(n) - \mathbf{w}^*(n) \rVert^2\) の推移をジャンプ時刻の前後で並べたもので、\(\lambda\) が小さいほどジャンプ直後の誤差が急速に減衰する(追従が速い)一方、ジャンプ前の定常誤差の底(ノイズによる揺らぎの幅)が持ち上がる様子が一目で分かります。

計算量とLMS/NLMSとの比較

https://yuhi-sa.github.io/posts/20260310_adaptive_filter/1/ で示した比較表を再掲します。

特性LMSNLMSRLS
計算量(1ステップ)\(O(M)\)\(O(M)\)\(O(M^2)\)
メモリ\(O(M)\)\(O(M)\)\(O(M^2)\)
収束速度遅い(線形)中程度速い(二次収束)
定常誤差ステップサイズ依存ステップサイズ依存小さい(\(\lambda\) 依存)
数値安定性高い高い\(P\) 行列の発散に注意

RLSは行列 \(P\) の更新を反復するため、有限精度演算では丸め誤差が蓄積して発散することがあります。この数値的な落とし穴と、それ以外の実務上の注意点を次節でまとめて扱います。

実務上の注意点:数値的な落とし穴とエッジケース

行列反転補題によるRLSの再帰更新は数学的には厳密ですが、実装・運用の場面では以下のエッジケースに注意が必要です。

1. \(P\) 行列の数値的発散。 有限精度演算では丸め誤差が蓄積し、理論上は対称正定値であるはずの \(P(n)\) が非対称・非正定値になって発散することがあります。対策として、毎ステップ \(P \leftarrow (P + P^\top)/2\) で対称化する、あるいはQR分解やUD分解(Bierman分解)に基づくRLS変種を使う方法が定番です。本記事のPython実装(式(5)をそのまま実装した素朴な形)は教育目的であり、長時間運用する組み込み用途ではこの対策が必須です。

2. 忘却係数の「ワインドアップ」問題。 \(\lambda < 1\) の場合、入力が一時的に乏しい区間(無入力・定数入力など)が続くと、新しい情報がないまま式(5)の \(1/\lambda\) 倍だけが繰り返され、\(P(n)\) が指数的に増大し続けます。これは「covariance windup(共分散の暴走)」と呼ばれる現象で、その後に強い入力が再開すると、肥大化した \(P\) に起因する過大なカルマンゲイン \(\mathbf{k}(n)\) のせいで係数が一時的に大きく暴れることがあります。対策としては、\(\operatorname{tr} P(n)\) に上限を設ける、方向性忘却(directional forgetting:励起されている方向にのみ忘却を適用する)を用いる、あるいは残差(イノベーション)の大きさに応じて \(\lambda\) を動的に調整する可変忘却係数RLS(variable forgetting factor RLS)に切り替える、といった手法が使われます。

3. 持続的励起条件(persistent excitation)の欠如。 入力ベクトル \(\mathbf{x}(n)\) がある部分空間内でしか変動しない(例:定数入力、単一正弦波、ランク落ちした回帰変数)場合、その方向については \(\mathbf{R}(n)\) が特異または悪条件になり、対応する係数成分は一意に決まりません。\(\lambda = 1\) ならその方向の推定値は単に凍結されるだけですが、\(\lambda < 1\) では新しい情報が入らないままその方向の不確実性(\(P\) の対応する固有値)だけが増大し続け、上記のワインドアップと同じ問題を引き起こします。設計上、フィルタの入力が全タップ方向を十分にカバーするように励起されているか(周波数成分が偏っていないか)を事前に確認すべきです。

4. 初期値 \(\delta\) の選び方。 \(P(0) = I / \delta\) の \(\delta\) は、\(n < M\) の初期段階で \(\mathbf{R}(n)\) がまだランク落ちしている(観測数がタップ数に満たない)ことに対する正則化項です。\(\delta\) を大きくしすぎる(\(P(0)\) を小さくしすぎる)と初期の係数推定がゼロに強く引っ張られてバイアスが生じ、逆に \(\delta\) を小さくしすぎる(\(P(0)\) を大きくしすぎる)と初期の数ステップで \(\mathbf{k}(n)\) が過大になり推定が振動します。本記事の実験のように入力を単位分散に正規化してある場合は \(\delta = 1\) 程度が扱いやすく、\(\delta\) の影響は \(n \gg M\) では実質的に消えます。

5. \(\lambda\) の実務的な範囲。 定常誤差と追従性のトレードオフ(前節参照)を踏まえると、実務では \(\lambda \in [0.95, 0.999]\) 程度の範囲で調整されることが多く、\(\lambda < 0.9\) まで下げるケースは非常に速い非定常性を追う必要がある特殊な用途に限られます。\(\lambda\) を下げすぎると、有効に使われるサンプル数(実効メモリ長は概ね \(1/(1-\lambda)\) )が過小になり、ノイズに対する感度が高くなりすぎる点に注意してください。

用途別ガイド

用途推奨理由
タップ数が小さく(\(M \lesssim 20\) )高速収束が必須RLS二次収束で最短時間で正規方程式の解に到達
タップ数が大きい・組み込み低リソースLMS / NLMS\(O(M)\) で計算・メモリコストを抑制
環境が緩やかに非定常RLS(\(\lambda < 1\) )忘却係数で定常誤差と追従性を調整可能
環境の変化が急激・大きいカルマンフィルタプロセスノイズ \(Q\) を明示的にモデル化できる
長時間連続運用(数値安定性が重要)NLMS\(P\) 行列の発散リスクがなく実装も単純

関連記事

おすすめ書籍

はじめて学ぶディジタル・フィルタと高速フーリエ変換(三上直樹、CQ出版)

適応フィルタを含むディジタルフィルタの基礎を体系的に学べる定番書です。本記事のRLS導出の理論的背景を補完します。

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

参考文献

  • Haykin, S. (2014). Adaptive Filter Theory (5th ed.). Pearson. Chapters 9-13.
  • Sayed, A. H. (2008). Adaptive Filters. Wiley-IEEE Press. Chapter 5 (RLS).
  • Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35-45.
  • Haykin, S. (2001). Kalman Filtering and Neural Networks. Wiley. Chapter 1 (RLS-Kalman equivalence).
  • Lai, B., & Bernstein, D. S. (2024). Adaptive Kalman Filtering Developed from Recursive Least Squares Forgetting Algorithms. arXiv:2404.10914.