再帰的最小二乗法

## 再帰的最小二乗法

最小二乗法の解は次のように書けた。notationは[こちらの記事](https://qiita.com/ground0state/items/74eb6b4beef6bd520f5c)に従う。

```math
\hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top \boldsymbol{y}
```

これはデータの計画行列で表されるため、行列のデータ数の次元を落とす操作により、$n$個あるデータのうち$k$個までを使用した推定値$\hat{\boldsymbol{\beta}}(k)$と、$k-1$個まで使用した推定値$\hat{\boldsymbol{\beta}}(k-1)$の漸化式を与えることが出来る。このようにして再帰的にパラメータを更新して最小二乗法の解を推定する方法を再帰的最小二乗法と呼ぶ。

$P\equiv (X^\top X)^{-1}$と定義する。$k$個までのデータで構成した計画行列$X(1:k)$を使用して、$k$個までのデータを含む計画行列で構成した$P$を$P(k)$と書くことにする。

```math
\begin{align}
\left( P^{-1}(k) \right)_{ij} 
&= \left( X^\top(1:k) X(1:k) \right)_{ij}\\
&=\sum_{m=1}^k X^\top_{im}(1:k) X_{mj}(1:k) \\
&= \sum_{m=1}^k x_{mi} x_{mj} \\
&= \sum_{m=1}^{k-1} x_{mi} x_{mj} + x_{ki} x_{kj} \\
&= \sum_{m=1}^{k-1} X^\top_{im} X_{mj} + X^\top_{ik} X_{kj} \\
&= \left( X^\top(1:k-1) X(1:k-1) \right)_{ij} + X^\top_{ik} X_{kj} \\
\therefore P(k) &= \left( P(k-1)^{-1} + \boldsymbol{x}(k) \boldsymbol{x}^\top(k) \right)^{-1}
\end{align}
```
ここで$\boldsymbol{x}(k)$は$k$番目のデータの特徴量ベクトルで$(p+1)$行$1$列ベクトルである。これが$P$の更新則である。

次に$X^\top \boldsymbol{y}$について考える。

```math
\begin{align}
X^\top(1:k) \boldsymbol{y}(1:k) &= \sum_{m=1}^{k} \boldsymbol{x}(m)y_m \\
&= \sum_{m=1}^{k-1} \boldsymbol{x}(m)y_m +  \boldsymbol{x}(k)y_k \\
&= X^\top(1:k-1) \boldsymbol{y}(1:k-1) +  \boldsymbol{x}(k)y_k \\
&= P(k-1)^{-1}\hat{\boldsymbol{\beta}}(k-1) +  \boldsymbol{x}(k)y_k
\end{align}
```
$P$の更新則と組み合わせて次の式を得る。

```math
\begin{align}
\hat{\boldsymbol{\beta}}(k) 
&= P(k)X^\top(1:k) \boldsymbol{y}(1:k) \\
&= P(k)\left( P(k-1)^{-1}\hat{\boldsymbol{\beta}}(k-1) +  \boldsymbol{x}(k)y_k \right) \\
&= P(k)\left( \left( P^{-1}(k) - \boldsymbol{x}(k)\boldsymbol{x}^\top(k) \right) \hat{\boldsymbol{\beta}}(k-1) +  \boldsymbol{x}(k)y_k \right) \\
&= P(k)\left( P^{-1}(k)\hat{\boldsymbol{\beta}}(k-1) - \boldsymbol{x}(k)\boldsymbol{x}^\top (k) \hat{\boldsymbol{\beta}}(k-1)  +  \boldsymbol{x}(k)y_k \right) \\
&=  \hat{\boldsymbol{\beta}}(k-1) - P(k)\boldsymbol{x}(k)\boldsymbol{x}^\top(k) \hat{\boldsymbol{\beta}}(k-1)  +  P(k)\boldsymbol{x}(k)y_k \\
&=  \hat{\boldsymbol{\beta}}(k-1) +  P(k)\boldsymbol{x}(k)\left(y_k - \boldsymbol{x}^\top(k) \hat{\boldsymbol{\beta}}(k-1) \right) \\
&=  \hat{\boldsymbol{\beta}}(k-1) +  e(k)P(k)\boldsymbol{x}(k) \\
\end{align}
```
ここで$e(k)$は$k-1$番目のデータで推定した$\hat{\boldsymbol{\beta}}(k-1)$を用いて、$k$番目のデータ$(\boldsymbol{x}(k), y_k)$を回帰予測した残差である。

```math
e(k) = y_k - \boldsymbol{x}^\top(k) \hat{\boldsymbol{\beta}}(k-1)
```

これで$ \hat{\boldsymbol{\beta}}$の更新式を得られた。

さらに、$P$の更新式に逆行列計算が含まれているので、逆行列を計算しない更新式に変形する。[Woodburyの恒等式](https://manabitimes.jp/math/1041)、

```math
(A+B D C)^{-1}=A^{-1}-A^{-1} B\left(D^{-1}+C A^{-1} B\right)^{-1} C A^{-1}
```
を使用して$A=P(k-1)^{-1}$、$B=\boldsymbol{x}(k)$、$C=\boldsymbol{x}^\top(k)$、$D=1$とすると、$P$の更新則は次のように書ける。

```math
\begin{align}
P(k) &= \left( P(k-1)^{-1} + \boldsymbol{x}(k) \boldsymbol{x}^\top(k) \right)^{-1} \\
&= P(k-1) - P(k-1)\boldsymbol{x}(k)\left(1 +  \boldsymbol{x}^\top(k)P(k-1)\boldsymbol{x}(k) \right)^{-1} \boldsymbol{x}^\top(k)P(k-1) \\
&= P(k-1) - \frac{P(k-1)\boldsymbol{x}(k)\boldsymbol{x}^\top(k)P(k-1)}{1 +  \boldsymbol{x}^\top(k)P(k-1)\boldsymbol{x}(k)}
\end{align}
```

以上より、$P$と$\hat{\boldsymbol{\beta}}$の更新式を得た。

ここで、$k=1$のときが問題になる。

```math
\begin{align}
\hat{\boldsymbol{\beta}}(1) &= \hat{\boldsymbol{\beta}}(0) +  e(1)P(1)\boldsymbol{x}(1) \\
P(1) &= P(0) - \frac{P(0)\boldsymbol{x}(1)\boldsymbol{x}^\top(1)P(0)}{1 +  \boldsymbol{x}^\top(1)P(0)\boldsymbol{x}(1)} \\
\end{align}
```

もし、今考えている$n$個のデータ以前にもデータが存在している場合は、それらを初期値にすることが出来る。

```math
\begin{align}
\hat{\boldsymbol{\beta}}(0) &= \beta_{\mathrm{old}} \\
P(0) &= (X^\top_{\mathrm{old}} X_{\mathrm{old}})^{-1}
\end{align}
```

存在しない場合はデータが$(p+1)$個まで溜まってから最小二乗法を行って初期値にするか、次のように設定することが多い。※要出典

```math
\begin{align}
\hat{\boldsymbol{\beta}}(0) &= \boldsymbol{0} \\
P(0) &= \alpha I
\end{align}
```

ここで$\alpha$は$10^3\sim 10^5$である。

一般には忘却係数を導入して、古いデータの重みを小さくすることもできる。

- [Recursive least squares filter](https://en.wikipedia.org/wiki/Recursive_least_squares_filter)
- [逐次最小二乗法（RLS: Recursive Least Squares）ーパラメータ同定に用いられる時系列データの回帰手法！ー](https://yuyumoyuyu.com/2020/12/27/recursiveleastsquares/)