会議用のスピーカーフォンには、スピーカーから出た自分側の音声が回り込んでマイクに入る「音響エコー」を消す回路が入っています。このエコー経路は、部屋の形・人の位置・机の上の書類の有無で刻々と変わります。だからフィルタ係数を工場で固定しておくことはできず、通話が始まった瞬間から数百ミリ秒で経路を推定し、誰かが立ち上がれば即座に追従し直す必要があります。ここで問題になるのが「どれだけ速く正解にたどり着けるか」です。
適応フィルタの定番であるLMSアルゴリズムは、誤差の瞬時値を勾配の代わりに使って係数をじわじわ動かします。実装は驚くほど軽い一方で、入力信号に強い相関がある（低域に偏っている、といった）場合には収束が極端に遅くなります。音声のように相関の強い信号では、LMSが正解に近づく前に環境のほうが変わってしまうこともあります。RLS（Recursive Least Squares、再帰最小二乗法） は、この「速さ」を正面から取りに行くアルゴリズムです。毎サンプル、これまでに観測した全データに対する最小二乗解を厳密に解き直す — それを $O(M^2)$ の再帰計算で実現します。
適応フィルタの仕事は、未知システムの中身を一切知らないまま、入力と出力（望ましい応答）だけを見てその入出力関係を真似ることです。図の破線が示すとおり、誤差 $e(n)$ が係数の修正量に変換されて適応フィルタへ戻る、というフィードバックループになっています。LMSとRLSの違いはこの「修正量の作り方」だけであり、構成図そのものは両者で完全に共通です。
RLSは応用範囲が広く、たとえば次の場面で使われます。
- 通信路の適応等化 — 高速移動体との無線リンクでは、短いトレーニング系列の間にチャネル逆特性を掴み切る必要があります。RLSはフィルタ長 $M$ の2倍程度のサンプル数で収束するため、トレーニング系列を大幅に短縮できます。
- 音響エコーキャンセラ・能動騒音制御 — 相関の強い音声・騒音信号でも収束が鈍らないため、立ち上がりの速さが体感品質に直結する用途で選ばれます。
- システム同定と故障診断 — プラントの伝達特性をオンライン推定し、忘却係数で「古いデータを忘れる」ことで、機器の劣化や急変を検出できます。
本記事の内容
- LMSとRLSの発想の違い（瞬時勾配 vs 各時刻での厳密最小化）
- 指数重み付き最小二乗規範 $J(n)$ と正規方程式の導出
- 相関行列の再帰式と、行列反転補題（Sherman–Morrison）による逆行列の再帰更新
- ゲインベクトル $\bm{k}(n)$・事前誤差 $\xi(n)$・更新式 $\bm{w}(n) = \bm{w}(n-1) + \bm{k}(n)\xi^*(n)$ の1行ずつの導出
- 初期化 $\bm{P}(0) = \delta^{-1}\bm{I}$ が持つ正則化の意味
- 忘却係数 $\lambda$ と実効メモリ $1/(1-\lambda)$、追従性と定常誤差のトレードオフ
- 計算量 $O(M^2)$ と数値安定性（$\lambda < 1$ での発散）
- Pythonによるシステム同定・適応等化の実測比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 適応フィルタとLMSアルゴリズム — 最急降下法による係数更新
- ウィーナフィルタの理論 — 最小平均二乗誤差の最適フィルタを導出する
- 等化器（イコライザ）の設計 — 符号間干渉を除去する受信信号処理
- ディジタルフィルタの基礎（FIR/IIR）を解説
LMSとRLSは何が違うのか
同じ「適応フィルタ」でも、LMSとRLSでは目指しているものが根本的に違います。この違いを最初に押さえておくと、以降の長い導出が何のためのものか見失わずに済みます。
霧の中で山を下ることを想像してください。LMSは、足元の傾きだけを頼りに一歩ずつ下る登山者です。一歩の幅（ステップサイズ $\mu$）を大きくすれば速く下れますが、大きすぎると谷を飛び越えて振動します。しかも足元の傾きは「その瞬間の1サンプルから推定した勾配」なので、ノイズだらけです。つまりLMSは、ノイズの乗った勾配を頼りに、小刻みに、遠回りをしながら谷底へ向かいます。
一方RLSは、これまで歩いてきた道の地形図をすべて持っている登山者です。1歩進むごとに地形図を更新し、「今持っている全データから見て、いま最も低い点はどこか」を毎回計算し直して、そこへ直接ワープします。だから最初の数歩で谷底の近くに着いてしまいます。代償は、地形図（正確には逆相関行列 $\bm{P}(n)$）を持ち歩き、毎ステップ更新するためのメモリと計算量です。
タップ数を2に減らして係数平面上に描くと、この比喩が文字どおりの絵になります。入力に相関があると誤差曲面の等高線は細長い楕円になり、LMS（赤）は谷の幅が広い方向にはすぐ動けるのに、谷に沿った細長い方向にはほとんど進めず、260サンプル経っても出発点の近くで足踏みしています。一方RLS（青）は数ステップで大きくジャンプし、谷底の星印（真の係数）にほぼ乗ってしまいます。この「等高線の細長さ」を数値化したものが後で出てくる固有値広がりです。
この対比を数式で書くと、両者の違いは「何を最小化しているか」に集約されます。LMSが暗黙に目指しているのは、統計的な期待値としての平均二乗誤差
$$ J_{\text{MSE}}(\bm{w}) = E\!\left[|d(n) – \bm{w}^H \bm{u}(n)|^2\right] $$
であり、その最小化に必要な勾配 $\nabla J = -2(\bm{p} – \bm{R}\bm{w})$ を、期待値を取らずに1サンプルで置き換えたのがLMSです。RLSが最小化するのは、期待値ではなく実際に観測されたデータの二乗和です。
$$ \begin{equation} J(n) = \sum_{i=1}^{n} \lambda^{n-i} \left| d(i) – \bm{w}^H(n)\, \bm{u}(i) \right|^2 \end{equation} $$
ここで $\bm{u}(i) = [u(i), u(i-1), \dots, u(i-M+1)]^T$ はタップ入力ベクトル、$d(i)$ は望ましい応答（教師信号）、$\bm{w}(n)$ は時刻 $n$ における長さ $M$ の係数ベクトルです。$(\cdot)^H$ はエルミート転置（共役転置）を表し、実信号なら単なる転置と読み替えて構いません。
見落としてはいけないのが、$\bm{w}(n)$ が総和の外にある点です。時刻 $n$ で決める1本の係数ベクトルを、過去の全時刻 $i = 1, \dots, n$ の誤差に共通して適用し、その重み付き二乗和を最小にします。つまりRLSは各時刻ごとに「これまでの全データを説明する最良のフィルタ」を求め直しているわけです。LMSのように前の係数から少し動かすのではなく、原理的には毎回ゼロから解き直しています。これが速さの正体です。
そして $\lambda^{n-i}$ という重みが付いています。$0 < \lambda \le 1$ の定数 $\lambda$ を忘却係数（forgetting factor）と呼びます。最新のデータ $i = n$ には重み $\lambda^0 = 1$、1つ前には $\lambda$、$k$ サンプル前には $\lambda^k$ が掛かるので、古いデータほど指数的に軽くなります。これは「昔の部屋の音響特性は、もう今の役に立たない」という時変環境への構えです。$\lambda = 1$ とすれば全データが等しく効き、普通の（重みなし）最小二乗法に戻ります。
では、この $J(n)$ を最小にする $\bm{w}(n)$ を実際に求めてみましょう。
指数重み付き最小二乗の正規方程式
まず $J(n)$ を展開します。誤差を $e(i) = d(i) – \bm{w}^H\bm{u}(i)$ と書くと、$|e(i)|^2 = e(i)e^*(i)$ なので
$$ |e(i)|^2 = \left(d(i) – \bm{w}^H\bm{u}(i)\right)\left(d^*(i) – \bm{u}^H(i)\bm{w}\right) $$
となります。$\bm{w}^H\bm{u}(i)$ はスカラーなので、その共役は $\bm{u}^H(i)\bm{w}$ です。この積を展開すると4つの項が出ます。
$$ |e(i)|^2 = |d(i)|^2 – \bm{w}^H \bm{u}(i) d^*(i) – d(i)\,\bm{u}^H(i)\bm{w} + \bm{w}^H \bm{u}(i)\bm{u}^H(i)\bm{w} $$
これに $\lambda^{n-i}$ を掛けて $i$ について足し上げます。$\bm{w}$ は $i$ に依存しないので総和の外へ出せて、
$$ \begin{equation} J(n) = r_d(n) – \bm{w}^H \bm{z}(n) – \bm{z}^H(n)\bm{w} + \bm{w}^H \bm{\Phi}(n)\bm{w} \end{equation} $$
という二次形式にまとまります。ここで登場した3つの量が、RLSの主役です。
$$ \begin{equation} \bm{\Phi}(n) = \sum_{i=1}^{n} \lambda^{n-i}\, \bm{u}(i)\bm{u}^H(i) \end{equation} $$
$$ \begin{equation} \bm{z}(n) = \sum_{i=1}^{n} \lambda^{n-i}\, \bm{u}(i)\, d^*(i) \end{equation} $$
そして $r_d(n) = \sum_i \lambda^{n-i}|d(i)|^2$ は $\bm{w}$ を含まない定数項です。$\bm{\Phi}(n)$ は入力の時間平均相関行列（$M \times M$、エルミート、半正定値）、$\bm{z}(n)$ は入力と望ましい応答の時間平均相互相関ベクトル（$M$ 次元）にあたります。ウィーナフィルタで出てきた $\bm{R} = E[\bm{u}\bm{u}^H]$ と $\bm{p} = E[\bm{u}d^*]$ の、期待値を指数重み付き時間平均に置き換えたもの、と理解すると見通しが良くなります。
最小化のために $\bm{w}$ で微分します。複素ベクトルの微分では、$\bm{w}$ と $\bm{w}^*$ を独立変数とみなして $\bm{w}^*$ で偏微分するWirtinger微分を使うのが定石です。$\partial(\bm{w}^H\bm{z})/\partial\bm{w}^* = \bm{z}$、$\partial(\bm{z}^H\bm{w})/\partial\bm{w}^* = \bm{0}$、$\partial(\bm{w}^H\bm{\Phi}\bm{w})/\partial\bm{w}^* = \bm{\Phi}\bm{w}$ を使うと、
$$ \frac{\partial J(n)}{\partial \bm{w}^*} = -\bm{z}(n) + \bm{\Phi}(n)\bm{w} = \bm{0} $$
を得ます。整理すると、これが正規方程式（normal equation）です。
$$ \begin{equation} \bm{\Phi}(n)\, \bm{w}(n) = \bm{z}(n) \end{equation} $$
実信号の場合も、$J(n)$ を各係数 $w_k$ で偏微分してゼロと置けば同じ式に到達します。$\bm{\Phi}(n)$ が正則なら解は
$$ \bm{w}(n) = \bm{\Phi}^{-1}(n)\, \bm{z}(n) $$
です。二次形式の係数行列 $\bm{\Phi}(n)$ が半正定値なので、この停留点は確かに最小値を与えます。
さて、原理はこれで終わりですが、このままでは実用になりません。毎サンプル $M \times M$ 行列を逆行列にすると $O(M^3)$ かかります。$M = 256$ タップのエコーキャンセラなら、1サンプルあたり1600万回規模の演算です。48 kHzサンプリングでは到底追いつきません。ここから先が、RLSがRLSたるゆえんの工夫 — 前の時刻の逆行列を使い回して、次の逆行列を $O(M^2)$ で作る、という話になります。
相関行列の再帰式
まず、逆行列の前に $\bm{\Phi}(n)$ 自身の再帰式を作ります。定義式 (3) から、$i = n$ の項だけを分離してみましょう。
$$ \bm{\Phi}(n) = \sum_{i=1}^{n-1} \lambda^{n-i}\bm{u}(i)\bm{u}^H(i) + \bm{u}(n)\bm{u}^H(n) $$
第1項の指数を $\lambda^{n-i} = \lambda \cdot \lambda^{(n-1)-i}$ と分解すると、残った総和はちょうど $\bm{\Phi}(n-1)$ の定義そのものです。したがって
$$ \begin{equation} \bm{\Phi}(n) = \lambda\, \bm{\Phi}(n-1) + \bm{u}(n)\bm{u}^H(n) \end{equation} $$
まったく同様に、相互相関ベクトルについても
$$ \begin{equation} \bm{z}(n) = \lambda\, \bm{z}(n-1) + \bm{u}(n)\, d^*(n) \end{equation} $$
が成り立ちます。式 (6) の意味は明快です。「前の相関行列を $\lambda$ 倍して少し忘れ、新しい観測 $\bm{u}(n)$ の外積を1つ足す」。この更新は $O(M^2)$ で済みます。
そして重要なのは、加える項 $\bm{u}(n)\bm{u}^H(n)$ がランク1だという点です。
更新の中身を絵にすると、やっていることは「全体をひと回り薄くする（$\lambda$ 倍＝忘れる）」と「新しい外積を1枚重ねる」の2手だけだと分かります。中央のヒートマップは $\bm{u}(n)$ の要素の積でできているため、行も列も互いに比例関係にあり、右端の棒グラフが示すとおり固有値は1個（4.29）を残して残りは厳密にゼロです。この「1方向の情報しか加わらない」という強い構造こそが、逆行列を安く更新できる根拠になります。ランク1の摂動を受けた行列の逆行列は、元の逆行列から陽に計算できます。それを保証するのが、次に見る行列反転補題です。
行列反転補題（Sherman–Morrison）
行列反転補題は、次の主張です。
補題: $\bm{A}$（$M\times M$）と $\bm{C}$（$N\times N$）が正則で、$\bm{B}$ が $M\times N$、$\bm{D}$ が $N\times M$ のとき、
$$ \begin{equation} (\bm{A} + \bm{B}\bm{C}\bm{D})^{-1} = \bm{A}^{-1} – \bm{A}^{-1}\bm{B}\left(\bm{C}^{-1} + \bm{D}\bm{A}^{-1}\bm{B}\right)^{-1}\bm{D}\bm{A}^{-1} \end{equation} $$
$N = 1$（$\bm{B}$ が列ベクトル、$\bm{D}$ が行ベクトル、$\bm{C}$ がスカラー）の場合が特にSherman–Morrisonの公式と呼ばれます。
証明は、右辺が本当に逆行列になっていることを直接掛け算で確かめるのが最も明快です。右辺を $\bm{X}$ と書き、$(\bm{A} + \bm{B}\bm{C}\bm{D})\bm{X} = \bm{I}$ を示します。表記を短くするため $\bm{S} = \bm{C}^{-1} + \bm{D}\bm{A}^{-1}\bm{B}$ と置きます。まず素直に展開すると、
$$ (\bm{A} + \bm{B}\bm{C}\bm{D})\bm{X} = \bm{I} – \bm{B}\bm{S}^{-1}\bm{D}\bm{A}^{-1} + \bm{B}\bm{C}\bm{D}\bm{A}^{-1} – \bm{B}\bm{C}\bm{D}\bm{A}^{-1}\bm{B}\bm{S}^{-1}\bm{D}\bm{A}^{-1} $$
の4項が出ます。第2〜4項をまとめるため、共通する $\bm{B}$ を左に、$\bm{D}\bm{A}^{-1}$ を右にくくり出すと、
$$ (\bm{A} + \bm{B}\bm{C}\bm{D})\bm{X} = \bm{I} + \bm{B}\left[-\bm{S}^{-1} + \bm{C} – \bm{C}\,\bm{D}\bm{A}^{-1}\bm{B}\,\bm{S}^{-1}\right]\bm{D}\bm{A}^{-1} $$
角括弧の中を見ます。$\bm{C}$ を $\bm{C}\bm{S}\bm{S}^{-1}$ と書き直せば、$\bm{C}\bm{S} = \bm{C}(\bm{C}^{-1} + \bm{D}\bm{A}^{-1}\bm{B}) = \bm{I} + \bm{C}\bm{D}\bm{A}^{-1}\bm{B}$ なので、
$$ -\bm{S}^{-1} + \left(\bm{I} + \bm{C}\bm{D}\bm{A}^{-1}\bm{B}\right)\bm{S}^{-1} – \bm{C}\bm{D}\bm{A}^{-1}\bm{B}\bm{S}^{-1} = \bm{0} $$
と3項がきれいに打ち消し合います。よって $(\bm{A} + \bm{B}\bm{C}\bm{D})\bm{X} = \bm{I}$ が示され、$\bm{X}$ は確かに逆行列です。$\square$
補題のご利益は、左辺の $M\times M$ の逆行列を、右辺では $N\times N$（ここでは $1\times 1$、つまりただの割り算）の逆行列で置き換えられる点にあります。ランク1更新なら、逆行列の更新はスカラー除算と行列ベクトル積だけで済むわけです。
RLSへの適用
式 (6) に補題を当てはめます。対応は
$$ \bm{A} = \lambda\bm{\Phi}(n-1), \quad \bm{B} = \bm{u}(n), \quad \bm{C} = 1, \quad \bm{D} = \bm{u}^H(n) $$
です。以後、逆相関行列を $\bm{P}(n) = \bm{\Phi}^{-1}(n)$ と書きます。$\bm{A}^{-1} = \lambda^{-1}\bm{P}(n-1)$ に注意して補題を適用すると、
$$ \begin{equation} \bm{P}(n) = \lambda^{-1}\bm{P}(n-1) – \frac{\lambda^{-2}\,\bm{P}(n-1)\bm{u}(n)\bm{u}^H(n)\bm{P}(n-1)}{1 + \lambda^{-1}\bm{u}^H(n)\bm{P}(n-1)\bm{u}(n)} \end{equation} $$
を得ます。分母は $1\times1$ のスカラーです。ここで、分数の中に何度も現れる $\bm{P}(n-1)\bm{u}(n)$ をひとまとめにした量を定義します。
$$ \begin{equation} \bm{k}(n) = \frac{\lambda^{-1}\bm{P}(n-1)\bm{u}(n)}{1 + \lambda^{-1}\bm{u}^H(n)\bm{P}(n-1)\bm{u}(n)} = \frac{\bm{P}(n-1)\bm{u}(n)}{\lambda + \bm{u}^H(n)\bm{P}(n-1)\bm{u}(n)} \end{equation} $$
2つ目の等号は、分子分母に $\lambda$ を掛けただけです。実装では右側の形（分母が $\lambda + \text{スカラー}$）を使うのが素直です。この $\bm{k}(n)$ をゲインベクトル（gain vector）と呼びます。定義を使うと式 (9) は
$$ \begin{equation} \bm{P}(n) = \lambda^{-1}\left[\bm{P}(n-1) – \bm{k}(n)\bm{u}^H(n)\bm{P}(n-1)\right] \end{equation} $$
という簡潔なリッカチ方程式になります。行列ベクトル積が2回程度、外積が1回で更新できるので、計算量は $O(M^2)$ です。
ゲインベクトルのもうひとつの顔
ここで、RLSの導出で最も気持ちのよい補題を証明します。
主張: $\bm{k}(n) = \bm{P}(n)\bm{u}(n)$。つまりゲインベクトルは、更新後の逆相関行列を現在の入力に掛けたものに等しい。
証明します。式 (11) の両辺に右から $\bm{u}(n)$ を掛けると、
$$ \bm{P}(n)\bm{u}(n) = \lambda^{-1}\bm{P}(n-1)\bm{u}(n) – \lambda^{-1}\bm{k}(n)\,\bm{u}^H(n)\bm{P}(n-1)\bm{u}(n) $$
です。ここで第1項を書き換えるために、式 (10) の定義を分母を払った形
$$ \lambda^{-1}\bm{P}(n-1)\bm{u}(n) = \bm{k}(n)\left[1 + \lambda^{-1}\bm{u}^H(n)\bm{P}(n-1)\bm{u}(n)\right] $$
に直して代入します。すると
$$ \bm{P}(n)\bm{u}(n) = \bm{k}(n) + \lambda^{-1}\bm{k}(n)\bm{u}^H(n)\bm{P}(n-1)\bm{u}(n) – \lambda^{-1}\bm{k}(n)\bm{u}^H(n)\bm{P}(n-1)\bm{u}(n) = \bm{k}(n) $$
となり、第2項と第3項が完全に打ち消し合って主張が示されます。$\square$
この事実が効いてくるのは次節です。$\bm{k}(n) = \bm{P}(n)\bm{u}(n) = \bm{\Phi}^{-1}(n)\bm{u}(n)$ という表現は、ゲインベクトルが「入力ベクトルを、これまでの入力の相関構造で白色化したもの」だと教えてくれます。入力にエネルギーが集中している方向（$\bm{\Phi}$ の大きい固有値の方向）は $\bm{\Phi}^{-1}$ で押さえ込まれ、逆にほとんど励振されていない方向は増幅されます。LMSが $\mu\bm{u}(n)$ という「生の入力方向」に係数を動かすのに対し、RLSは相関で正規化した方向に動かす — これが固有値広がりの影響を受けない理由の芽です。
タップ数2で実際に描くと、この「白色化」が目に見えます。$\rho = 0.95$ の入力は右上がりの帯状に分布しており、これが $\bm{\Phi}$ の大きい固有値の方向です。同じ入力 $\bm{u}(n)$ に対し、LMSは入力そのままの $17^\circ$ 方向へ動くのに対し、RLSのゲインは $-42^\circ$ — ほぼ帯と直交する方向を向きます。何度も観測して十分に情報のある方向は縮められ、めったに励振されない方向が持ち上げられているわけです。
RLSの係数更新式を導く
準備が整いました。ゴールは、$\bm{w}(n) = \bm{P}(n)\bm{z}(n)$ を「$\bm{w}(n-1)$ からの補正」の形に書き直すことです。
出発点は正規方程式の解と、相互相関ベクトルの再帰式 (7) です。
$$ \bm{w}(n) = \bm{P}(n)\bm{z}(n) = \bm{P}(n)\left[\lambda\,\bm{z}(n-1) + \bm{u}(n)d^*(n)\right] $$
括弧を展開して2つの項に分けます。
$$ \bm{w}(n) = \lambda\,\bm{P}(n)\bm{z}(n-1) + \bm{P}(n)\bm{u}(n)\,d^*(n) $$
第2項には、いま証明した $\bm{P}(n)\bm{u}(n) = \bm{k}(n)$ がそのまま使えて、$\bm{k}(n)d^*(n)$ になります。第1項は、リッカチ方程式 (11) の両辺に $\lambda$ を掛けた $\lambda\bm{P}(n) = \bm{P}(n-1) – \bm{k}(n)\bm{u}^H(n)\bm{P}(n-1)$ を代入します。
$$ \lambda\bm{P}(n)\bm{z}(n-1) = \bm{P}(n-1)\bm{z}(n-1) – \bm{k}(n)\bm{u}^H(n)\bm{P}(n-1)\bm{z}(n-1) $$
ここで $\bm{P}(n-1)\bm{z}(n-1) = \bm{w}(n-1)$（1時刻前の正規方程式の解）ですから、右辺は $\bm{w}(n-1) – \bm{k}(n)\bm{u}^H(n)\bm{w}(n-1)$ です。2つの項を合わせると、
$$ \bm{w}(n) = \bm{w}(n-1) – \bm{k}(n)\bm{u}^H(n)\bm{w}(n-1) + \bm{k}(n)d^*(n) $$
$\bm{k}(n)$ でくくれば、
$$ \bm{w}(n) = \bm{w}(n-1) + \bm{k}(n)\left[d^*(n) – \bm{u}^H(n)\bm{w}(n-1)\right] $$
角括弧の中身に名前を付けましょう。事前誤差（a priori estimation error）を
$$ \begin{equation} \xi(n) = d(n) – \bm{w}^H(n-1)\bm{u}(n) \end{equation} $$
と定義すると、その共役は $\xi^*(n) = d^*(n) – \bm{u}^H(n)\bm{w}(n-1)$ で、角括弧そのものです。したがって最終形は
$$ \begin{equation} \bm{w}(n) = \bm{w}(n-1) + \bm{k}(n)\,\xi^*(n) \end{equation} $$
となります。「事前誤差」という名前は、まだ更新していない古い係数 $\bm{w}(n-1)$ で予測したときの誤差だからです。実装上も、$\bm{w}$ を更新する前に計算できる量なので都合が良いのです。
式 (13) の形はLMSの更新式 $\bm{w}(n) = \bm{w}(n-1) + \mu\,\bm{u}(n)e^*(n)$ とそっくりです。違いはただ一点、方向ベクトルが $\mu\bm{u}(n)$（生の入力にスカラーのステップサイズ）から $\bm{k}(n) = \bm{\Phi}^{-1}(n)\bm{u}(n)$（相関で白色化された入力）に変わっただけ。RLSは「入力相関に合わせてステップサイズを方向ごとに自動調整したLMS」だと見ることもできます。しかもその調整はヒューリスティックではなく、最小二乗解に厳密に一致するように決まっています。
事後誤差と変換係数
もうひとつ、RLSの安定性を直感的に支える性質を見ておきます。更新後の係数で測った誤差を事後誤差（a posteriori error）$e(n) = d(n) – \bm{w}^H(n)\bm{u}(n)$ と呼びます。式 (13) のエルミート転置 $\bm{w}^H(n) = \bm{w}^H(n-1) + \xi(n)\bm{k}^H(n)$ を代入すると、
$$ e(n) = d(n) – \bm{w}^H(n-1)\bm{u}(n) – \xi(n)\bm{k}^H(n)\bm{u}(n) = \xi(n)\left[1 – \bm{k}^H(n)\bm{u}(n)\right] $$
を得ます。$\bm{k}(n) = \bm{P}(n)\bm{u}(n)$ かつ $\bm{P}(n)$ はエルミートなので $\bm{k}^H(n)\bm{u}(n) = \bm{u}^H(n)\bm{P}(n)\bm{u}(n)$ です。そこで
$$ \begin{equation} \gamma(n) = 1 – \bm{u}^H(n)\bm{P}(n)\bm{u}(n) = \frac{\lambda}{\lambda + \bm{u}^H(n)\bm{P}(n-1)\bm{u}(n)} \end{equation} $$
と置くと $e(n) = \gamma(n)\,\xi(n)$ です（2つ目の等号は、$\bm{u}^H\bm{P}(n)\bm{u} = \bm{u}^H\bm{k}(n)$ に式 (10) を代入すれば出ます）。この $\gamma(n)$ を変換係数と呼びます。
$\bm{P}(n-1)$ が正定値なら分母は $\lambda$ より大きいので、$0 < \gamma(n) \le 1$ が常に成り立ちます。つまり事後誤差は必ず事前誤差より小さい。更新は必ず誤差を減らす方向にしか動かない、ということです。$\gamma(n)$ はまた「その入力方向をどれだけ既に学習済みか」の指標でもあります。まったく新しい方向の入力が来ると $\bm{u}^H\bm{P}\bm{u}$ が大きくなり $\gamma \to 0$、すなわち事後誤差ほぼゼロまで一気に合わせ込みます。すでに何度も見た方向なら $\gamma \approx 1$ で、ほとんど動きません。
実測でも $\gamma(n)$ は最小 $0.006$・最大 $0.976$ と、理論どおり $0$ と $1$ の間に完全に収まりました。左のグラフを見ると、学習し始めの十数サンプルは $\gamma \approx 0$ で「新しい方向にはためらわず一気に合わせ込む」挙動を示し、その後は $\gamma$ が $1$ に張り付いて更新がほとんど止まります。右のグラフでは事後誤差（青）が全サンプルで事前誤差（赤）以下に収まっており（実測で100.0%）、水色に塗った差分が1回の更新で減らせた誤差にあたります。学習初期はこの差が2桁近くありますが、収束後はほぼゼロです。
以上で、RLSの全構成要素が揃いました。次はこれを実行可能な手順にまとめ、まだ触れていない初期値の問題を片付けます。
アルゴリズムの全体像と初期化
RLSの1サンプルあたりの手順は次の5行です。
- ゲイン計算: $\displaystyle \bm{k}(n) = \frac{\bm{P}(n-1)\bm{u}(n)}{\lambda + \bm{u}^H(n)\bm{P}(n-1)\bm{u}(n)}$
- 事前誤差: $\xi(n) = d(n) – \bm{w}^H(n-1)\bm{u}(n)$
- 係数更新: $\bm{w}(n) = \bm{w}(n-1) + \bm{k}(n)\xi^*(n)$
- 逆相関行列更新: $\bm{P}(n) = \lambda^{-1}\left[\bm{P}(n-1) – \bm{k}(n)\bm{u}^H(n)\bm{P}(n-1)\right]$
- 対称化（実装上の保険）: $\bm{P}(n) \leftarrow \frac{1}{2}\left[\bm{P}(n) + \bm{P}^H(n)\right]$
手順1と4で $\bm{P}(n-1)\bm{u}(n)$ を共有できるので、行列ベクトル積は実質1回で済みます。
$\bm{P}(0) = \delta^{-1}\bm{I}$ は何をしているのか
さて、$n = 0$ では観測がまだ1つもないので $\bm{\Phi}(0) = \bm{0}$、その逆行列は存在しません。困ります。実務では
$$ \begin{equation} \bm{P}(0) = \delta^{-1}\bm{I}, \qquad \bm{w}(0) = \bm{0} \end{equation} $$
と初期化します（$\delta > 0$ は小さな正の定数）。これは単なる場当たり的な逃げではなく、はっきりした意味があります。$\bm{P}(0) = \delta^{-1}\bm{I}$ は $\bm{\Phi}(0) = \delta\bm{I}$ を意味しますから、再帰式 (6) を $n$ 回適用すると
$$ \bm{\Phi}(n) = \sum_{i=1}^{n}\lambda^{n-i}\bm{u}(i)\bm{u}^H(i) + \delta\lambda^{n}\bm{I} $$
になります。つまりRLSが実際に解いているのは、規範
$$ J(n) = \sum_{i=1}^{n}\lambda^{n-i}|e(i)|^2 + \delta\lambda^{n}\|\bm{w}(n)\|^2 $$
の最小化 — リッジ回帰（$L_2$正則化つき最小二乗） です。データがまだ少なく $\bm{\Phi}$ が特異になりがちな初期に、係数が暴れるのを正則化項が押さえます。しかも正則化の強さ $\delta\lambda^n$ は時間とともに指数的に減衰するので、データが溜まれば自動的に正則化が抜けて、素の最小二乗解に近づきます。よくできた設計です。
$\delta$ の選び方には定石があります。$\delta$ を小さく（$\bm{P}(0)$ を大きく）すれば初期のゲイン $\bm{k}$ が大きくなり、立ち上がりが速くなりますが、雑音の多い環境では初期に暴れます。逆に $\delta$ を大きくすると滑らかですが立ち上がりが鈍ります。目安として、高SNRなら $\delta$ を入力分散の $0.01$ 倍程度、低SNRなら $1$ 倍程度に取るのが一つの指針です。
雑音を多めにした条件（40回平均）で $\delta$ を4段階に振ると、この定石が実測で確認できます。$n=20$ でのミスアライメントは $\delta = 0.001$ で $-18.8$ dB、$\delta = 0.01$ で $-18.9$ dB、$\delta = 1$ で $-14.7$ dB、$\delta = 100$ では $-4.7$ dB。$\delta$ を小さくするほど立ち上がりが速いことがはっきり出ています。ただし $\delta = 0.001$（赤）は最初の数サンプルで一度 $-1$ dB 付近まで持ち上がっており、これが「正則化が弱すぎて初期に暴れる」現象です。一方 $\delta = 100$（濃紺）は正則化が強すぎて $200$ 反復後でも $-11.6$ dB にとどまり、他の3本が到達している $-29$ dB にはるかに及びません。
初期化まで含めてアルゴリズムが確定したので、次は残るハイパーパラメータ $\lambda$ の意味を掘り下げます。
忘却係数 $\lambda$ と実効メモリ
$\lambda$ の役割は「どれだけ過去を覚えているか」の一点に尽きます。重み $\lambda^{n-i}$ の総和を取ると、$n \to \infty$ で
$$ \sum_{k=0}^{\infty}\lambda^{k} = \frac{1}{1-\lambda} $$
です。この値は「重み1のデータが何個分あるのに相当するか」を表すので、実効メモリ長（有効窓長）と呼ばれます。$\lambda = 0.99$ なら100サンプル、$\lambda = 0.95$ なら20サンプル、$\lambda = 0.999$ なら1000サンプル分です。データの「平均年齢」で測るなら $\sum_k k\lambda^k(1-\lambda) = \lambda/(1-\lambda)$ で、実効メモリ長よりちょうど1だけ小さい値になります。
左のグラフは、$\lambda = 0.95$ なら100サンプル前のデータの重みがすでに $0.006$ とほぼ無視できる大きさまで落ちるのに対し、$\lambda = 0.999$ では400サンプル前でもまだ $0.67$ の重みが残ることを示しています。右のグラフの縦軸が対数目盛であることに注目してください。$\lambda$ を $0.95 \to 0.99 \to 0.999$ とわずかに $1$ へ寄せるだけで、記憶するデータ数は $20 \to 100 \to 1000$ と1桁ずつ跳ね上がります。$\lambda$ の設計が「小数第3位まで気を遣う」作業になるのはこのためです。
この一本の数字が、RLSの2大性能を正反対の向きに支配します。
定常誤差（ミスアジャストメント）: $\lambda$ が小さいほど、実効的に使うデータ数が減るので、雑音による推定のばらつきが大きくなります。過剰平均二乗誤差を最小平均二乗誤差 $J_{\min}$ で割ったミスアジャストメントは、よく知られた近似として
$$ \begin{equation} \mathcal{M} \approx M\,\frac{1 – \lambda}{1 + \lambda} \end{equation} $$
と表されます（$M$ はフィルタ長）。$\lambda \to 1$ でミスアジャストメントはゼロに近づき、理論上は最適解に完全収束します。
追従性: 一方、$\lambda \to 1$ にすると過去を全部覚えてしまうので、環境が変わったとき、古いデータの重みに引きずられて追従が遅れます。$\lambda$ を小さくすれば古い情報を素早く捨て、新しい環境にすぐ乗り換えられます。
つまり $\lambda$ は「精度 vs 俊敏さ」のつまみです。定常環境なら $\lambda = 0.999$ 以上、時変環境なら $0.95 \sim 0.99$、という選び方になります。実効メモリ長を「環境がほぼ一定とみなせる時間」に合わせる、と考えるのが実践的な設計指針です。
計算量のトレードオフ
RLSの代償ははっきりしています。
| LMS | NLMS | RLS | |
|---|---|---|---|
| 1サンプルの乗算回数 | $\approx 2M$ | $\approx 3M$ | $\approx 4M^2$ |
| メモリ | $O(M)$ | $O(M)$ | $O(M^2)$ |
| 収束（反復数の目安） | 固有値広がりに強く依存 | 依存が緩和 | $\approx 2M$、広がりに非依存 |
| 数値安定性 | 高い | 高い | 注意が必要 |
左の両対数グラフでは、LMS/NLMSが傾き1の直線であるのに対しRLSだけが傾き2の直線になり、$M$ が伸びるほど差が開いていくことが視覚的に分かります。$M = 8$ ではLMSの16倍で済みますが、$M = 256$ では512倍です。右は実測時間で、行列反転補題を使うRLS更新（青）と、毎サンプル $\bm{\Phi}(n)$ を作って np.linalg.solve で解き直す素朴な実装（紫）の比較です。$M = 256$ で $182$ マイクロ秒 対 $1549$ マイクロ秒と約9倍の差がつきました。$O(M^3)$ 対 $O(M^2)$ の理論差ほど開かないのは、BLASの行列演算が高度に最適化されているためですが、それでも補題を使う価値は明確です。
$M = 256$ なら、LMSの512回に対しRLSは約26万回。500倍以上の差です。この差があるため、長いフィルタが必要な音響エコーキャンセラでは、素のRLSではなく、$\bm{P}$ の構造（入力ベクトルのシフト構造）を利用して $O(M)$ に落とした高速RLS（FTF、格子型RLSなど）が使われます。ただし高速化した版は数値的にさらに脆く、周期的なリセットや安定化項が必須になります。
計算量の話が出たところで、RLSのもうひとつの弱点 — 数値安定性 — を見ておきましょう。
数値安定性: $\lambda < 1$ での発散
RLSの $\bm{P}(n)$ は、理論上は常に正定値エルミート行列です。ところが有限語長の演算では、リッカチ方程式 (11) の引き算で誤差が蓄積し、$\bm{P}$ の対称性が崩れたり、固有値が負になったりします。しかも式 (11) には $\lambda^{-1}$ という 1より大きい倍率が毎ステップ掛かります。$\lambda = 0.95$ なら1ステップで $1/0.95 \approx 1.053$ 倍。誤差が減衰せず、指数的に増幅されていくわけです。
特に危険なのが励振不足（persistent excitation の欠如）の状況です。入力が止まる、あるいは入力が単一正弦波でフィルタ長より少ない次元しか張らない場合、$\bm{\Phi}(n)$ のある固有方向にはもう新しい情報が加わりません。式 (6) を見れば明らかで、そんな方向では $\bm{\Phi}$ の固有値が毎ステップ $\lambda$ 倍に縮み続け、$\bm{P}$ の対応する固有値は $\lambda^{-1}$ 倍に膨らみ続けます。これが有名な共分散ワインドアップ（covariance wind-up、$\bm{P}$ の爆発）です。$\bm{P}$ が天文学的な値になった状態でわずかな雑音が入ると、係数が一瞬で吹き飛びます。
対策はいくつかあります。
- $\lambda = 1$ にする（時不変なら） — $\lambda^{-1}$ の増幅がなくなり、$\bm{P}$ は単調に小さくなるので発散しません。ただし追従性は失われます。
- 可変忘却係数 — 誤差が小さいとき（環境が定常なとき）は $\lambda \to 1$、誤差が大きくなったら $\lambda$ を下げる。
- 正則化リーク — $\bm{P}(n) \leftarrow \bm{P}(n) + \epsilon\bm{I}$ の代わりに $\bm{\Phi}(n) \leftarrow \bm{\Phi}(n) + \epsilon\bm{I}$ 相当の処理を入れ、$\bm{P}$ の固有値に上限を設ける。$\text{tr}(\bm{P})$ が閾値を超えたらリセットする実装もよく使われます。
- 平方根版RLS — $\bm{P}$ を直接持たず、そのコレスキー因子 $\bm{P} = \bm{L}\bm{L}^H$ を更新する（QR分解RLS、逆QR-RLS）。因子で持てば正定値性が構造的に保証され、条件数も平方根に改善します。実機では第一選択です。
理屈はここまでです。実際にどれだけ違うのか、数値で確かめましょう。
Pythonでの実装
RLSとLMS/NLMSの実装
まずアルゴリズム本体を実装します。RLSは前節の5行の手順をそのまま書き下すだけです。比較対象としてLMSとNLMSも用意します。
import numpy as np
def tap_vector(x, n, M):
"""時刻nでのタップ入力ベクトル u(n) = [x(n), x(n-1), ..., x(n-M+1)]"""
u = x[n::-1][:M]
if len(u) < M: # 先頭では0埋め
u = np.concatenate([u, np.zeros(M - len(u))])
return u
def run_rls(x, d, M, lam=0.99, delta=0.01):
"""RLS適応フィルタ（実数版）。係数履歴と事前誤差を返す"""
N = len(x)
P = np.eye(M) / delta # P(0) = delta^-1 I（リッジ正則化）
w = np.zeros(M) # w(0) = 0
W = np.zeros((N, M))
err = np.zeros(N)
for n in range(N):
u = tap_vector(x, n, M)
Pu = P @ u # 1回だけ計算して使い回す
k = Pu / (lam + u @ Pu) # ゲインベクトル 式(10)
xi = d[n] - w @ u # 事前誤差 式(12)
w = w + k * xi # 係数更新 式(13)
P = (P - np.outer(k, u @ P)) / lam # リッカチ方程式 式(11)
P = 0.5 * (P + P.T) # 対称化（数値安定化）
W[n], err[n] = w, xi
return W, err
RLSの本体はループ内のわずか5行です。行列演算は P @ u（$O(M^2)$）と np.outer(k, u @ P)（$O(M^2)$）だけで、逆行列計算はどこにも現れません。これが行列反転補題のご利益です。続いて比較用のLMS/NLMSを実装します。
import numpy as np
def run_lms(x, d, M, mu=0.01, normalized=False, eps=1e-6):
"""LMS / NLMS 適応フィルタ"""
N = len(x)
w = np.zeros(M)
W = np.zeros((N, M))
err = np.zeros(N)
for n in range(N):
u = tap_vector(x, n, M)
e = d[n] - w @ u # 誤差
if normalized: # NLMS: 入力パワーで正規化
w = w + mu * e * u / (eps + u @ u)
else: # LMS: 生の入力方向へ
w = w + mu * e * u
W[n], err[n] = w, e
return W, err
LMSの更新はスカラー倍と加算だけで $O(M)$、RLSと比べると圧倒的に軽いことがコードの見た目からも分かります。次に、この差が収束速度としてどう現れるかを実験で確かめます。
実験1: システム同定と固有値広がりの影響
未知のFIRシステム $\bm{w}_o$（長さ $M = 8$）に既知の入力を通し、その出力に微小な観測雑音を加えたものを $d(n)$ とします。適応フィルタが $\bm{w}_o$ をどれだけ速く突き止められるかを、係数誤差ノルム $\|\bm{w}(n) – \bm{w}_o\|^2$（ミスアライメント）で測ります。
入力には AR(1) 過程 $x(n) = \rho\,x(n-1) + \sqrt{1-\rho^2}\,v(n)$ を使います。$\rho$ を上げるほど入力の相関が強くなり、相関行列 $\bm{R}$ の固有値広がり $\chi = \lambda_{\max}/\lambda_{\min}$ が大きくなります。ここがLMSの弱点を突く設定です。
import numpy as np
def make_ar1(N, rho, rng):
"""分散1に正規化したAR(1)入力"""
v = rng.standard_normal(N)
x = np.zeros(N)
x[0] = v[0]
for n in range(1, N):
x[n] = rho * x[n - 1] + np.sqrt(1 - rho ** 2) * v[n]
return x
def eigen_spread(M, rho):
"""AR(1)入力の相関行列 R[i,j] = rho^|i-j| の固有値広がり"""
i = np.arange(M)
R = rho ** np.abs(i[:, None] - i[None, :])
ev = np.linalg.eigvalsh(R)
return ev.max() / ev.min()
for rho in [0.0, 0.9, 0.95]:
print(f"rho={rho}: 固有値広がり chi = {eigen_spread(8, rho):.1f}")
出力は rho=0.0: chi = 1.0、rho=0.9: chi = 113.4、rho=0.95: chi = 263.8 です。白色入力では相関行列が単位行列なので広がりは1、$\rho = 0.95$ では260倍以上に達します。LMSの収束時定数は最小固有値に支配されるため、この広がりがそのまま「最も遅い方向の遅さ」になります。RLSはゲインに $\bm{\Phi}^{-1}$ が入っているので、この広がりが打ち消されるはずです。実際に走らせて確かめます。
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
def ensemble(rho, N=2000, M=8, trials=100, sigma=0.01, seed=0):
"""アンサンブル平均のミスアライメント曲線を返す"""
wo = np.random.default_rng(seed).standard_normal(M) / np.sqrt(M)
curves = {"LMS (μ=0.01)": [], "NLMS (μ=0.2)": [],
"RLS (λ=0.99)": [], "RLS (λ=0.95)": []}
for t in range(trials):
r = np.random.default_rng(1000 + t)
x = make_ar1(N, rho, r)
U = np.stack([tap_vector(x, n, M) for n in range(N)])
d = U @ wo + sigma * r.standard_normal(N)
Ws = [run_lms(x, d, M, 0.01)[0], run_lms(x, d, M, 0.2, normalized=True)[0],
run_rls(x, d, M, 0.99, 0.01)[0], run_rls(x, d, M, 0.95, 0.01)[0]]
for key, W in zip(curves, Ws):
curves[key].append(np.sum((W - wo) ** 2, axis=1))
return {k: np.mean(v, axis=0) for k, v in curves.items()}
100回の独立試行で平均を取っているのは、単一試行の学習曲線は雑音で暴れて比較にならないためです。適応フィルタの性能評価では、このアンサンブル平均が標準的な手続きです。次に結果を描画します。
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(16, 4.5), sharey=True)
for ax, rho in zip(axes, [0.0, 0.9, 0.95]):
for name, m in ensemble(rho).items():
ax.plot(10 * np.log10(m), linewidth=1.4, label=name)
ax.axhline(-20, color="gray", linestyle=":", linewidth=1)
ax.set_title(f"入力相関 ρ={rho}（固有値広がり χ={eigen_spread(8, rho):.0f}）")
ax.set_xlabel("反復回数 n")
ax.grid(alpha=0.3)
axes[0].set_ylabel("係数誤差 $\\|w(n)-w_o\\|^2$ [dB]")
axes[0].legend(fontsize=9)
plt.tight_layout()
plt.show()
3枚を並べると、RLSとLMSの性格の違いが一目で分かります。RLSの2本（濃青系）は $\rho$ を変えても左端でほぼ垂直に落ちる形が変わりません。対してLMS（赤）は右の図に行くほど傾きが寝ていき、$\rho = 0.95$ では2000反復かけても $-20$ dB の点線に届いていません。定常状態の高さ（グラフ右端の水平部分）を見ると、収束の速いRLS($\lambda$=0.95)がRLS($\lambda$=0.99)より一段高い位置で止まっており、速さと精度が引き換えであることも読み取れます。実測値を整理すると次の通りです。
| 入力 | 固有値広がり $\chi$ | LMSが $-20$ dBに到達 | NLMS | RLS($\lambda$=0.99) | RLS($\lambda$=0.95) |
|---|---|---|---|---|---|
| 白色 ($\rho$=0) | 1.0 | 201 反復 | 85 | 10 | 10 |
| AR(1) ($\rho$=0.9) | 113 | 1445 反復 | 239 | 10 | 10 |
| AR(1) ($\rho$=0.95) | 264 | 到達せず | 358 | 10 | 10 |
この表を対数目盛の棒グラフにすると、RLSの2本だけが3条件を通じて高さ10で完全に横並びなのに対し、LMSとNLMSの棒は右へ行くほど伸びていくのが際立ちます。斜線の棒は2000反復以内に到達しなかったことを示す打ち切り表示です。読み取れることは3つあります。第一に、RLSの収束はフィルタ長 $M = 8$ に対して10反復 — つまり理論が言う「$2M$ 反復程度」とほぼ一致し、しかも入力相関がどれだけ強くなってもまったく変化しません。ゲインベクトルに $\bm{\Phi}^{-1}(n)$ が入っているため、入力の固有値広がりが構造的に打ち消されているのです。第二に、LMSは白色入力の201反復から $\rho = 0.9$ で1445反復へと7倍以上に悪化し、$\rho = 0.95$ では2000反復以内に $-20$ dBへ届きません。第三に、NLMSは入力パワーによる正規化のおかげでLMSよりずっと健闘しますが（85 → 239 → 358反復）、それでもRLSとは1桁以上の差が残ります。パワーによる正規化はスカラーの調整に過ぎず、方向ごとの調整はできないからです。
実験2: 忘却係数と定常誤差のトレードオフ
同じ実験の定常状態に注目すると、$\lambda$ の効果が見えます。$\rho = 0.9$ での最終200サンプルの平均ミスアライメントは、$\lambda = 0.99$ で $-44.5$ dB、$\lambda = 0.95$ で $-37.2$ dB でした。差は $7.3$ dB です。
これは理論と合っているでしょうか。式 (16) のミスアジャストメントは $\lambda$ にしか依らない比を与えるので、2つの $\lambda$ の比を取ると
$$ \frac{\mathcal{M}(0.95)}{\mathcal{M}(0.99)} = \frac{(1-0.95)/(1+0.95)}{(1-0.99)/(1+0.99)} = \frac{0.02564}{0.005025} = 5.10 $$
すなわち $10\log_{10}5.10 = 7.08$ dB の差が予測されます。実測 $7.3$ dB とよく一致しています。$\lambda$ を小さくして忘却を強めた分だけ、実効的なデータ数が減り、雑音による推定のばらつきが増えているわけです。
では、$\lambda$ を大きくしておけば常に得なのでしょうか。それを確かめるのが次の実験です。
実験3: プラント急変時の追従性
$n = 1000$ で未知システムを別のものに差し替え、それぞれのアルゴリズムがどれだけ速く再収束するかを見ます。
import numpy as np
import matplotlib.pyplot as plt
def ensemble_switch(rho=0.9, N=2000, M=8, trials=100, sigma=0.01, change_at=1000):
"""n=change_at で未知システムが切り替わる場合の追従性"""
rg = np.random.default_rng(0)
wo1 = rg.standard_normal(M) / np.sqrt(M)
wo2 = rg.standard_normal(M) / np.sqrt(M)
Wtrue = np.tile(wo1, (N, 1))
Wtrue[change_at:] = wo2
curves = {"LMS (μ=0.01)": [], "NLMS (μ=0.2)": [],
"RLS (λ=0.99)": [], "RLS (λ=0.95)": []}
for t in range(trials):
r = np.random.default_rng(1000 + t)
x = make_ar1(N, rho, r)
U = np.stack([tap_vector(x, n, M) for n in range(N)])
d = np.einsum("ij,ij->i", U, Wtrue) + sigma * r.standard_normal(N)
Ws = [run_lms(x, d, M, 0.01)[0], run_lms(x, d, M, 0.2, normalized=True)[0],
run_rls(x, d, M, 0.99, 0.01)[0], run_rls(x, d, M, 0.95, 0.01)[0]]
for key, W in zip(curves, Ws):
curves[key].append(np.sum((W - Wtrue) ** 2, axis=1))
return {k: np.mean(v, axis=0) for k, v in curves.items()}
plt.figure(figsize=(9, 5))
res = ensemble_switch()
for name, m in res.items():
plt.plot(10 * np.log10(m), linewidth=1.4, label=name)
seg = 10 * np.log10(m[1000:])
n20 = np.argmax(seg < -20) if (seg < -20).any() else -1
print(f"{name:14s} 急変後に-20dBへ再到達: {n20} サンプル")
plt.axvline(1000, color="k", linestyle="--", linewidth=1)
plt.text(1010, 5, "プラント急変", fontsize=10)
plt.xlabel("反復回数 n"); plt.ylabel("係数誤差 [dB]")
plt.title("プラント急変時の追従性（入力相関 ρ=0.9）")
plt.legend(fontsize=9); plt.grid(alpha=0.3)
plt.tight_layout(); plt.show()
左のグラフの $n=1000$ で全手法の誤差が一斉に $0$ dB 付近まで跳ね上がっているのが、システムが別物に差し替わった瞬間です。注目すべきはその後の下り方で、切り替え前は最も低い位置にいたRLS($\lambda$=0.99)（濃い青緑）が、切り替え後はRLS($\lambda$=0.95)（濃紺）に大きく後れを取っています。右の横棒グラフがその差を数字にしたもので、出力は「RLS($\lambda$=0.95): 76サンプル、RLS($\lambda$=0.99): 304サンプル、NLMS: 476サンプル、LMS: 再到達せず」でした。実験2とは順位が逆転しています。定常精度では $\lambda = 0.99$ が $7.3$ dB 優位でしたが、急変への追従では $\lambda = 0.95$ が4倍速く復帰します。$\lambda = 0.95$ の実効メモリ長は20サンプルなので、切り替わりから20サンプルも経てば古いシステムの情報はほぼ消えています。一方 $\lambda = 0.99$ は100サンプル分を覚えているため、しばらく古い解に引きずられるわけです。
これが $\lambda$ の設計そのものです。「どのくらいの周期で環境が変わるか」を見積もり、実効メモリ長 $1/(1-\lambda)$ をその時定数に合わせる。定常精度が欲しければ長く、俊敏さが欲しければ短く。どちらも欲しい場合は、誤差が大きい間だけ $\lambda$ を下げる可変忘却係数を使います。
実験4: 励振不足による $\bm{P}$ の発散
最後に、RLSの数値的な弱点を実演します。最初の500サンプルだけ入力を与え、その後は入力を止めます。式 (6) から、入力が来なければ $\bm{\Phi}(n) = \lambda\bm{\Phi}(n-1)$ なので $\bm{P}$ は毎ステップ $\lambda^{-1}$ 倍に膨らむはずです。
import numpy as np
M, delta, N = 8, 0.01, 3000
for lam in [1.0, 0.99, 0.95]:
r = np.random.default_rng(3)
x = np.zeros(N)
x[:500] = r.standard_normal(500) # 前半のみ励振、後半は無入力
P, trace_log = np.eye(M) / delta, []
for n in range(N):
u = x[n::-1][:M]
if len(u) < M:
u = np.concatenate([u, np.zeros(M - len(u))])
Pu = P @ u
k = Pu / (lam + u @ Pu)
P = (P - np.outer(k, u @ P)) / lam
trace_log.append(np.trace(P))
print(f"λ={lam}: trace(P) n=500 → {trace_log[499]:.3e}, "
f"n=1500 → {trace_log[1499]:.3e}, n=2999 → {trace_log[-1]:.3e}")
縦軸が $10^{-2}$ から $10^{58}$ までの対数目盛であることに注意してください。入力が止まった $n=500$ を境に、$\lambda < 1$ の2本は対数目盛の上でまっすぐな直線になります。対数目盛での直線は指数増大を意味し、その傾きはちょうど $-\log\lambda$ です。結果は劇的です。$\lambda = 1.0$ では $\text{tr}(\bm{P})$ は $1.60\times10^{-2}$ のまま3000サンプルを通じて一定です。$\lambda = 0.99$ では $n=1500$ で $1.92\times10^{3}$、$n=2999$ で $6.77\times10^{9}$。$\lambda = 0.95$ に至っては $n=1500$ で $7.28\times10^{21}$、$n=2999$ で $1.89\times10^{55}$ という天文学的な値まで膨らみます。$\lambda^{-(n-500)}$ の指数増大がそのまま現れているわけです。
この状態でわずかでも雑音が乗った観測が入ると、巨大な $\bm{P}$ を通じてゲイン $\bm{k}$ が爆発し、係数は一瞬で意味を失います。実機で「しばらく無音が続いたあとエコーキャンセラが破綻する」といった不具合は、たいていこの共分散ワインドアップです。対策は前述の通りで、入力パワーを監視して励振が不足したら更新を止める、$\text{tr}(\bm{P})$ に上限を設けてリセットする、あるいは平方根版RLSを使う、といった実装上の防御が欠かせません。$\lambda = 1$ でも安全とは言え、そちらは追従性を捨てることになります。
RLSをいつ選ぶか
ここまでの理論と実測をまとめると、選択の指針は次のようになります。
RLSが向く場面: フィルタ長 $M$ が数十タップ以下で、収束の速さが決定的に重要な用途。短いトレーニング系列しか使えない通信路等化、瞬時にモデルを掴む必要があるオンライン同定、入力に強い相関がありLMSでは実用的な時間で収束しない場合。$O(M^2)$ の計算量が許容できるなら、RLSはほぼ常にLMS/NLMSより速く収束します。
LMS/NLMSが向く場面: フィルタ長が数百〜数千タップに及ぶ音響エコーキャンセラや能動騒音制御。演算資源が厳しい組み込み・固定小数点実装。数値的な堅牢性が最優先される長時間連続稼働。収束は遅くても、$O(M)$ の軽さと「壊れない」ことのほうが価値が高い場面は現実には多くあります。
中間の選択肢: アフィン射影法（APA）は、直近 $K$ 本の入力ベクトルだけを使って部分的に白色化するアルゴリズムで、$K=1$ でNLMS、$K=M$ でRLSに漸近します。$K$ を数本に取れば $O(KM)$ 程度の計算量で、NLMSとRLSの中間の収束速度が得られます。また、$\bm{P}$ をコレスキー因子で持つQR分解RLSは、計算量は $O(M^2)$ のまま数値安定性を大幅に改善します。
まとめ
本記事では、RLS（再帰最小二乗）適応フィルタを規範の定義から実装まで導出しました。
- 規範の違いがすべての出発点 — LMSが平均二乗誤差の勾配を1サンプルで推定するのに対し、RLSは指数重み付き二乗和 $J(n) = \sum_i \lambda^{n-i}|d(i) – \bm{w}^H\bm{u}(i)|^2$ を各時刻で厳密に最小化する
- 正規方程式 — $J(n)$ を展開して $\bm{w}^*$ で微分すると $\bm{\Phi}(n)\bm{w}(n) = \bm{z}(n)$ が得られ、$\bm{\Phi}(n) = \lambda\bm{\Phi}(n-1) + \bm{u}(n)\bm{u}^H(n)$ というランク1更新の再帰式が成り立つ
- 行列反転補題が計算量を落とす — ランク1更新にSherman–Morrisonを適用すると、$O(M^3)$ の逆行列計算が $O(M^2)$ のリッカチ方程式 $\bm{P}(n) = \lambda^{-1}[\bm{P}(n-1) – \bm{k}(n)\bm{u}^H(n)\bm{P}(n-1)]$ に化ける
- ゲインベクトルの正体 — $\bm{k}(n) = \bm{P}(n)\bm{u}(n) = \bm{\Phi}^{-1}(n)\bm{u}(n)$ は「相関で白色化された入力」であり、更新式 $\bm{w}(n) = \bm{w}(n-1) + \bm{k}(n)\xi^*(n)$ はLMSと同形でありながら固有値広がりの影響を受けない
- 初期化はリッジ正則化 — $\bm{P}(0) = \delta^{-1}\bm{I}$ は罰則項 $\delta\lambda^n\|\bm{w}\|^2$ を加えることと等価で、正則化の強さは時間とともに指数的に抜けていく
- 忘却係数はつまみ1本のトレードオフ — 実効メモリ長 $1/(1-\lambda)$ が長いほど定常誤差は小さく（実測で $\lambda$=0.99 は 0.95 より $7.3$ dB 良い）、短いほど追従が速い（急変後の再収束が4倍速い）
- 速さの代償 — 計算量 $O(M^2)$、メモリ $O(M^2)$、そして $\lambda<1$ では励振不足時に $\bm{P}$ が $\lambda^{-n}$ で発散する。実測では3000サンプルで $\text{tr}(\bm{P})$ が $10^{55}$ に達した
- 実測での収束差 — フィルタ長8に対しRLSは固有値広がり1でも264でも一貫して10反復で $-20$ dB に到達。同じ条件でLMSは201反復・1445反復・到達せず、と大きく劣化した
次のステップとして、以下の記事も参考にしてください。