26
20

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Last updated at Posted at 2016-09-04

スパースモデリング第3章追跡アルゴリズムのアルゴリズムをPythonで実装.

コードと実験結果をまとめたJupyter notebook

以下の貪欲法と凸緩和の手法を実装.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$に追加
  2. 暫定界の更新:$\| \mathbf{Ax}-\mathbf{b}\|^{2}_{2} , \rm{subject} , \rm{to} , \rm{Support}\{\mathbf{x}\}=S^{k}$の最適解を$\mathbf{x}^{k}$とする.
  3. 残差の更新:$\mathbf{r}^{k}=\mathbf{b}-\mathbf{Ax}^{k}$を計算する.
  4. 停止条件:もし$\|\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$を得る.
  2. 重みの更新:$X_{k}(j,j)=|x_{k}(j)|^{0.5}$として重み対角行列$\mathbf{X}$を更新
  3. 停止条件:もし$\|\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

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

平均サポート間距離
dist_S.png

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

X_1.png
X_3.png
X_5.png

考察

  • OMP結構頑張ってる
  • IRLSは教科書と結果が一致しないので原因の検討が必要
26
20
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
26
20

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?