スパースモデリング：第3章 追跡アルゴリズム

[スパースモデリング](http://www.kyoritsu-pub.co.jp/bookdetail/9784320123946)第3章追跡アルゴリズムのアルゴリズムをPythonで実装．

[コードと実験結果をまとめたJupyter notebook](https://github.com/kibo35/sparse-modeling/blob/master/ch03.ipynb)．

以下の貪欲法と凸緩和の手法を実装．IRLSはちょっと怪しい…

## 貪欲法
* 直交マッチング追跡(orthogonal matching pursuit; OMP)
* マッチング追跡(matching pursuit; MP)
* 弱マッチング追跡(weak matching pursuit; WMP)
* しきい値アルゴリズム

## 凸緩和の手法
* 反復再重み付け最小二乗法(iterative-reweighted-lease-sqaures; IRLS)

## 貪欲法の概要
#### 初期化
$k=0$として

* 初期解 $\mathbf{x}^{0}=\mathbf{0}$
* 初期残差 $\mathbf{r}^{0}=\mathbf{b}$
* 解の初期サポート$S^{0}=\emptyset$

#### メインループ 
$k \leftarrow k+1$として以下のステップを実行する．

1. 誤差の計算，サポートの更新：$\mathbf{A}$の中から一番残差$\mathbf{r}^{k}$を減らせる列$\mathbf{a}_{j}$を$S^k$に追加
1. 暫定界の更新：$\\| \mathbf{Ax}-\mathbf{b}\\|^{2}_{2} \, \rm{subject} \, \rm{to} \, \rm{Support}\\{\mathbf{x}\\}=S^{k}$の最適解を$\mathbf{x}^{k}$とする．
1. 残差の更新：$\mathbf{r}^{k}=\mathbf{b}-\mathbf{Ax}^{k}$を計算する．
1. 停止条件：もし$\\|\mathbf{r}^{k}\\|<\epsilon_{0}$なら終了し，そうでなければ反復する．

誤差の計算，サポートの更新では$\mathbf{A}$の全ての列$\mathbf{a}_{j}$と残差$\mathbf{r}$の内積を計算し，内積の絶対値が最大となる列をサポートに追加する．つまり$\mathbf{r}$に平行な列を順番に追加していく．

## 貪欲法のアルゴリズムの違い

* OMP
    * 上のまま
* MP
    * $\mathbf{A}$ではなく新たにサポートに追加した$\mathbf{a}_j$のみで暫定解を更新
* WMP
    * スカラー$t \, (0<t<1)$をハイパーパラメータに追加
    * 誤差の計算とサポートの更新で$\mathbf{r}^{k-1}$との内積の絶対値が$t\\|\mathbf{r}^{k-1}\\|$以上の$\mathbf{a}_j$が見つかったら計算を打ち切る．後はMPと同じ．
* しきい値アルゴリズム
    * サポートに含まれる基底関数の個数$k$をハイパーパラメータに追加
    * $\mathbf{b}$との内積の絶対値の大きい$\mathbf{a}_j$の上から$k$個をサポートとしてしまう．後は最小二乗問題を解き$\mathbf{x}$とする．


## IRLSの概要

* $\left( P_{0} \right)$問題ではなく$\left( P_{1} \right)$問題を解く．
* 重み付き$L_{2}$ノルムで$L_{1}$ノルムを近似して$\left( P_{1} \right)$問題を解く．
* スパースな解を得るためにはむっちゃ反復する必要がある．

#### 初期化
$k=0$として

* 初期近似解 $\mathbf{x}^{0}=\mathbf{1}$（全部1のベクトル）
* 初期重み行列 $\mathbf{X}^{0}=\mathbf{I}$

#### メインループ 
$k \leftarrow k+1$として以下のステップを実行する．

1. 最小二乗法：線形連立方程式$$\mathbf{x}\_{k} = \mathbf{X}\_{k-1}^{2} \mathbf{A}^{T} ( \mathbf{A} \mathbf{X}_{k-1}^{2} \mathbf{A}^{T} )^{+} \mathbf{b}$$を直接または反復法を用いて解き，近似解$\mathbf{x}\_k$を得る．
1. 重みの更新：$X_{k}(j,j)=\|x_{k}(j)\|^{0.5}$として重み対角行列$\mathbf{X}$を更新
1. 停止条件：もし$\\|\mathbf{x}\_{k}-\mathbf{x}\_{k-1}\\|\_{2}$が事前に与えられたしきい値より小さければ終了し，そうでなければ反復する．

実験
==

## シミュレーション
$\mathbf{x}$ 50次元のベクトル
$\mathbf{b}$ 30次元のベクトル
$\mathbf{A}$ 30×50次元の行列，$[-1, 1)$の一様乱数を列正規化

* $\mathbf{x}$の非ゼロの要素数$k$として、その非ゼロ要素の値を$[-2, -1] \cup [1, 2]$の範囲の一様分布からiid抽出する．
* $k=1,...,10$で各200回シミュレーション

## 指標

* 推定した$\mathbf{\hat{x}}$と真の$\mathbf{x}$の$L_2$ノルム
* 推定した$\mathbf{\hat{x}}$と真の$\mathbf{x}$のサポート間距離$$dist(\hat{S},S)=\frac{max\\{|\hat{S}|, |S| \\} - |\hat{S}\cap S|}{max\\{|\hat{S}|,|S|\\}}$$
* 200回のシミュレーションの平均をとる．

## 結果

平均$L_{2}$誤差
![l2_err.png](https://qiita-image-store.s3.amazonaws.com/0/114365/8cae6d8a-0831-f8c3-10e3-b73d64abdc89.png)

* kが大きくなると誤差が増える．
* IRLSの誤差が小さいけども$L_{2}$誤差だから当然か
* OMPも結構頑張ってる．速いし

平均サポート間距離
![dist_S.png](https://qiita-image-store.s3.amazonaws.com/0/114365/2241bfc7-8ffc-76f6-357b-a562b1c23d33.png)

* kが低いときのIRLSのサポート間距離が大きくなってしまった．教科書ではそんなことは無いので反復が足りないか、不具合かも
* OMPもがんばってる．

![X_1.png](https://qiita-image-store.s3.amazonaws.com/0/114365/6cd0cd44-a1a3-bd4b-04ad-aed4dc105804.png)
![X_3.png](https://qiita-image-store.s3.amazonaws.com/0/114365/a8c6e5d3-ab40-5eb5-2fc8-a978cda47fa5.png)
![X_5.png](https://qiita-image-store.s3.amazonaws.com/0/114365/e0d19af1-56d0-2994-348e-606c3afbf9eb.png)

考察
==
* OMP結構頑張ってる
* IRLSは教科書と結果が一致しないので原因の検討が必要