逐次最小二乗法(RLS)

逐次最小二乗法(Recursive Least Squares, RLS)について,問題設定から以下の更新式の導出まで解説します。

逐次最小二乗法における更新式

θ→n+1=θ→n+Pn+1a→n+1(bn+1−a→n+1⊤θ→n)Pn+1=Pn−Pna→n+1a→n+1⊤Pn1+a→n+1⊤Pna→n+1\begin{aligned} \overrightarrow{\theta}_{n+1} &= \overrightarrow {\theta}_{n}+P_{n+1}\overrightarrow{a}_{n+1} (b_{n+1}-\overrightarrow{a}_{n+1}^{\top}\overrightarrow{\theta}_n)\\ P_{n+1}&=P_n-\dfrac{P_n\overrightarrow{a}_{n+1}\overrightarrow{a}_{n+1}^{\top}P_n}{1+\overrightarrow{a}_{n+1}^{\top}P_n\overrightarrow{a}_{n+1}} \end{aligned}

逐次最小二乗法の問題設定

行列 AA と列ベクトル b→\overrightarrow{b} が与えられたときに,∥Aθ→−b→∥22\|A\overrightarrow{\theta}-\overrightarrow{b}\|_2^2 を最小にする θ→\overrightarrow{\theta} を求める問題を考えます(最小二乗法)。

この問題の重要性は,以下でも紹介しています。

AA は特徴量の値を並べた行列,b→\overrightarrow{b} は正解ラベル,θ→\overrightarrow{\theta} は推定したい dd 次元のパラメータです。

さらに,ここではデータが1つずつ追加されていく状況を考えます。つまり,

An=(a→1,a→2,…,a→n)⊤,b→n=(b1,b2,…,bn)⊤A_n=(\overrightarrow{a}_1,\overrightarrow{a}_2,\dots,\overrightarrow{a}_n)^{\top},\overrightarrow{b}_n=(b_1,b_2,\dots,b_n)^{\top}

という nn 個のデータがある場合の最小二乗パラメータ θn\theta_n は分かっている状況で,新しいデータ (a→n+1,bn+1)(\overrightarrow{a}_{n+1},b_{n+1}) を追加した場合のパラメータ θn+1\theta_{n+1} の計算方法を考えます。

最小二乗パラメータを計算する2つの方法

方法1:愚直に計算

n+1n+1 個のデータがある場合の最小二乗パラメータを愚直に直接計算すると,

θ→n+1=(An+1⊤An+1)−1An+1⊤b→n+1\overrightarrow{\theta}_{n+1}=(A_{n+1}^{\top}A_{n+1})^{-1}A_{n+1}^{\top}\overrightarrow{b}_{n+1}

となります。→最小二乗法の行列表現(単回帰,多変数,多項式)

方法2:逐次最小二乗法(RLS)

一方,1つ前の最適パラメータ θ→n\overrightarrow{\theta}_n と新しいデータ (a→n+1,bn+1)(\overrightarrow{a}_{n+1},b_{n+1}) から

θ→n+1=θ→n+Pn+1a→n+1(bn+1−a→n+1⊤θ→n)\overrightarrow{\theta}_{n+1}=\overrightarrow{\theta}_{n}+P_{n+1}\overrightarrow{a}_{n+1}(b_{n+1}-\overrightarrow{a}_{n+1}^{\top}\overrightarrow{\theta}_n)

という更新式で θ→n+1\overrightarrow{\theta}_{n+1} を計算することもできます。ただし,Pn+1P_{n+1} とは (An+1⊤An+1)−1(A_{n+1}^{\top}A_{n+1})^{-1} のことで,この d×dd\times d 行列も

Pn+1=Pn−Pna→n+1a→n+1⊤Pn1+a→n+1⊤Pna→n+1P_{n+1}=P_n-\dfrac{P_n\overrightarrow{a}_{n+1}\overrightarrow{a}_{n+1}^{\top}P_n}{1+\overrightarrow{a}_{n+1}^{\top}P_n\overrightarrow{a}_{n+1}}

という別の式を使って更新できます。つまり,θn\theta_n だけでなく PnP_n も覚えておく必要があります。

計算量の比較

方法1は,行列積や逆行列の計算を含むので計算に時間がかかりますが,方法2は O(d2)O(d^2) で計算できます。

このように「毎回1から全部計算し直す」のではなく「1つ前の結果をうまく更新して新しい結果を得る」ことで計算時間や使用メモリ量を削減できることがあります(オンラインアルゴリズム,逐次アルゴリズム)。

RLS の更新式の導出

θ→n+1=(An+1⊤An+1)−1An+1⊤b→n+1=Pn+1An+1⊤b→n+1\overrightarrow{\theta}_{n+1}=(A_{n+1}^{\top}A_{n+1})^{-1}A_{n+1}^{\top}\overrightarrow{b}_{n+1}\\=P_{n+1}A_{n+1}^{\top}\overrightarrow{b}_{n+1}

を,

θ→n=(An⊤An)−1An⊤b→n=PnAn⊤b→n\overrightarrow{\theta}_{n}=(A_{n}^{\top}A_{n})^{-1}A_{n}^{\top}\overrightarrow{b}_{n}=P_nA_{n}^{\top}\overrightarrow{b}_{n}

を使って表すのが目標です。逆行列の補助定理(Woodburyの恒等式)が登場するのが楽しいです。

導出

A⊤b→A^{\top}\overrightarrow{b} の部分を分解すると,

θ→n=Pn∑i=1na→ibi\overrightarrow{\theta}_n=P_n\displaystyle\sum_{i=1}^n\overrightarrow{a}_ib_i

となる。よって,

∑i=1na→ibi=Pn−1θ→n\displaystyle\sum_{i=1}^n\overrightarrow{a}_ib_i=P_n^{-1}\overrightarrow{\theta}_n

よって,

θ→n+1=Pn+1∑i=1n+1a→ibi=Pn+1(a→n+1bn+1+∑i=1na→ibi)=Pn+1(a→n+1bn+1+Pn−1θ→n)\overrightarrow{\theta}_{n+1}=P_{n+1}\displaystyle\sum_{i=1}^{n+1}\overrightarrow{a}_ib_i\\ =P_{n+1}(\overrightarrow{a}_{n+1}b_{n+1}+\displaystyle\sum_{i=1}^n\overrightarrow{a}_ib_i)\\ =P_{n+1}(\overrightarrow{a}_{n+1}b_{n+1}+P_n^{-1}\overrightarrow{\theta}_n)

ここで,

Pn+1−1=An+1⊤An+1=An⊤An+a→n+1a→n+1⊤=Pn−1+a→n+1a→n+1⊤P_{n+1}^{-1}=A_{n+1}^{\top}A_{n+1}\\ =A_{n}^{\top}A_n+\overrightarrow{a}_{n+1}\overrightarrow{a}_{n+1}^{\top}\\ =P_n^{-1}+\overrightarrow{a}_{n+1}\overrightarrow{a}_{n+1}^{\top}

という関係式(※)を使って上式の Pn−1P_n^{-1} を消去すると,

θ→n+1=Pn+1{a→n+1bn+1+(Pn+1−1−a→n+1a→n+1⊤)θ→n}=θ→n+Pn+1a→n+1(bn+1−a→n+1⊤θ→n)\overrightarrow{\theta}_{n+1}\\= P_{n+1}\{\overrightarrow{a}_{n+1}b_{n+1}+(P_{n+1}^{-1}-\overrightarrow{a}_{n+1}\overrightarrow{a}_{n+1}^{\top})\overrightarrow{\theta}_n\}\\ =\overrightarrow{\theta}_{n}+P_{n+1}\overrightarrow{a}_{n+1}(b_{n+1}-\overrightarrow{a}_{n+1}^{\top}\overrightarrow{\theta}_n)

となり,θ\theta の更新式を得た。

さらに,PnP_n の更新式を得るために,さきほどの関係式(※)に逆行列の補助定理

(A+BDC)−1=A−1−A−1B(D−1+CA−1B)−1CA−1(A+BDC)^{-1}\\=A^{-1}-A^{-1}B(D^{-1}+CA^{-1}B)^{-1}CA^{-1}

を使うと,

Pn+1=(Pn−1+a→n+1a→n+1⊤)−1=Pn−Pna→n+1(1+a→n+1⊤Pna→n+1)−1a→n+1⊤Pn=Pn−Pna→n+1a→n+1⊤Pn(1+a→n+1⊤Pna→n+1)P_{n+1}=(P_n^{-1}+\overrightarrow{a}_{n+1}\overrightarrow{a}_{n+1}^{\top})^{-1}\\ =P_n-P_n\overrightarrow{a}_{n+1}(1+\overrightarrow{a}_{n+1}^{\top}P_n\overrightarrow{a}_{n+1})^{-1}\overrightarrow{a}_{n+1}^{\top}P_n\\ =P_n-\dfrac{P_n\overrightarrow{a}_{n+1}\overrightarrow{a}_{n+1}^{\top}P_n}{(1+\overrightarrow{a}_{n+1}^{\top}P_n\overrightarrow{a}_{n+1})}

なお,A⊤AA^{\top}A に逆行列が存在しない場合のことは考えていません。実際,A⊤AA^{\top}A に逆行列が存在するくらい十分なデータがある状態からRLSを開始すればOKです。または,L2正則化(リッジ回帰)をすれば,A⊤A→A⊤A+λIA^{\top}A\to A^{\top}A+\lambda I となり,逆行列は必ず存在します。

全然高校数学ではないですが,逆行列の補助定理の応用例をぜひとも紹介したいので書きました。