スパースモデリングを超初歩から始めてみた。(ISTAを実装した)

ブログを書くのかなり久しぶりな感じがしますね。

春休み、ちょっと勉強しようと思って東北大の大関正之先生の書かれたスパースモデリングの入門テキスト

「今日からできるスパースモデリング」 大関正之

を読んで、すごいおもしろいなーと感じたので、ちょっとやってみることにしました。

次のような問題を考えます。\(A\)を\(M\times N\)行列として\(N\)次元ベクトル\(\boldsymbol{x}\)についての線形連立方程式
\[\boldsymbol{y}=A\boldsymbol{x}\]
を考えます。

このとき\(N < M\)であれば与えられた方程式の数の方が多いので方程式は解けます。これを優決定系といいます。

これに対して\(N>M\)であるときは条件式が足りないので方程式は解けないということになりますが、

ここで\(\boldsymbol{x}\)に次のような条件が事前に与えられているとしましょう。

\(K\)スパース性:\(\boldsymbol{x}\)のほとんどの成分は\(0\)であり,非ゼロ成分の個数は\(K\)個である。

\(K< M\)であれば方程式を解くことができます。

しかしながら現実の場面においては非ゼロ成分がどこに分布してるのかということや\(K\)自体がわからない場合が多いです。

このような場合、次のような最適化問題を解くことになりますが、これはNP困難な問題となり計算量爆発を起こしてしまいます。

\[\min_{\boldsymbol{x}}\|x\|_0\,\,\, \mathrm{s.t.}\,\, \boldsymbol{y}=A\boldsymbol{x}\]

ここで、\(\|\cdot\|_0\)は\(L_0\)ノルムと呼ばれ、非ゼロ成分の個数を表します。

そのような場合に有効なアプローチとして圧縮センシングの方法があるそうです。

これは\(L_0\)ノルムの代わりに\(L_1\)ノルム\(\|\cdot\|_1\)を用いるものです。\(L_1\)ノルムとは

\[\|\boldsymbol{x}\|=|x_1|+|x_2|+\ldots+|x_N|\]

で定義されるノルムのことです。

\(L_1\)ノルムを用いて次の最適化問題を考えます。

\(L_1\)再構成: \(\min_{\boldsymbol{x}}\|x\|_1\,\,\, \mathrm{s.t.}\,\, \boldsymbol{y}=A\boldsymbol{x}\)


\(L_1\)ノルムは非ゼロ成分の個数が少ないほど小さくなりますが、成分の大きさを反映しているため\(L_0\)ノルムによる最適化とは本質的には異なるはずですが、多くの場合うまくいくらしいです。

\(L_1\)再構成を具体的に解く計算アルゴリズムの基本的なものにISTA(Iterative Shrinkage Soft-thresholding Algorithm)というのがあります。次のようなアルゴリズムです。
ISTA

  1. \(t=0\),\(\boldsymbol{x} [0]\)を適当に初期化する。例えば\(\boldsymbol{x} [0]=A^\top \boldsymbol{y}\)。
  2. \(\boldsymbol{v}[t]=x[t]+(L\lambda)^{-1}A^\top(\boldsymbol{y}-A\boldsymbol{x}[t])\)
  3. \(\boldsymbol{x}[t+1]=S_{L^{-1}}(\boldsymbol{v}[t])\),\(t=t+1\)
  4. 終了条件を満たすまで2.ー3.を繰り返す。

ISTAの導出などについては上記の資料を見てください。\(S_\lambda(\cdot)\)というのは軟判定閾値関数とよばれるもので

\[
S_\lambda(y)=
\begin{cases}
y-\lambda &&y>\lambda\\
0 &&-\lambda\leq y\leq\lambda\\
y+\lambda &&y<-\lambda
\end{cases}
\]

で定義されます。\(S_\lambda(\boldsymbol{y})\)は\(\boldsymbol{y}\)の各成分に\(S_\lambda\)を作用させたものという意味です。

以下がコードです。なんだかあまり上手く原情報ベクトルが再現されてなかったのでどこかで間違ってそうです。

間違いなど、ドシドシ教えてください。それではまた。

import numpy as np
import random

#軟判定閾値関数
def ST(y,lam):
    if y>lam :
        return y-lam
    elif y<-lam :
        return y+lam
    else:
        return 0
N=1000
M=100
lam=1E-10
ep=1E-8
#原情報ベクトルの非零成分の個数
K=20
# 観測行列を生成
A=np.random.randn(M,N)
A_T=A.transpose()
#リプシッッツ定数の設定 L=|A^TA|_F(フロベニウスノベル)
L=np.linalg.norm(A_T.dot(A),ord="fro")/lam
print("L=%f" % L)
# スパースな原信号ベクトルを生成
x_0=np.zeros(N)
r_K=random.sample(range(N),K)
for i in r_K:
    x_0[i]=np.random.rand()
#観測ベクトルを生成
y=A.dot(x_0)
#ISTAの初期値の設定
t=0
x_t=A_T.dot(y)
#収束するまで繰り返す
while np.linalg.norm(y-A.dot(x_t)) > ep:
    v_t=x_t+A_T.dot(y-A.dot(x_t))/L/lam
    print("|y-Ax[%d]|_2=%f" % (t,np.linalg.norm(y-A.dot(x_t))))
    for i in range(N):
        x_t[i]=ST(v_t[i],1/L)
    t=t+1
print(np.linalg.norm(x_t-x_0))