スパースモデリングを超初歩から始めてみた。(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))

2015の振り返り(いまさら)

かなり久々の更新となります。
久しぶりすぎていったいどうしたものかという感じではあります。
更新してなかったのは書くべきネタがなかったからです。
本当になにもなかったのか、2015年のことをゆるゆる振り返りながら考えたいと思います。

2015年4月

希望に満ちて京都大学情報学科数理工学コースに配属されました。
数理工学コースのほうでは歓迎会があり、おいしい料理とともに数理の教員の方々に迎えていただきなんてすばらしいコースなんだと舌鼓を打っておりました。(食い気しかない)
計算機コースのほうでは歓迎会は開かれなかったようで、そこにもにんまりと優越感を感じていました。意地が悪いですね。
今年の数理の歓迎会は会場がカンフォーラ(京大正門を入って左手のレストラン)だったようで、さらに豪華な感じになってるんですかね?
いやあ、素晴らしい!おいしい料理に勝る歓迎はないと思います。*1

あとは必修科目の数理工学実験の登録を完全に忘れてて、担当教員に土下座メールかました苦い経験を思い出しました。
ほんとに後輩の皆さんは履修関係だけはほんとに気を付けるべきですよ!

2015年5月~7月

去年から継続してた朝永量子論ゼミをやったりしただけ。
特筆するべきことはなし。

2015年 8月~9月
夏休み。バイト漬けだった。

2015年 10月〜11月
サークルの定期演奏会に向けて練習してた。

2015年 12月 友達とHackU 2015というYahoo!主催のハッカソンに参加しました。メンバー全員がアプリ開発初心者だったので、完全にチャレンジ(笑)という状態でのスタートでした。本当になんとか発表前日までに完成に持って行けたのは留学生のKくんのおかげだと心から思います。Iくんの作ってくれたデモPVも素晴らしくて、そういう技を持ってる人はすごいなあという一言につきました。持っているというか、やる過程で頑張って勉強してくれたのかもしれませんが。どちらにせよすごいですね。僕はほとんど貢献できなかったので、悔しさが残るところですが、みんなでわいわいアイデアを形にしていく過程はとても楽しかったです。もう少しプログラミングの勉強をしとけば開発面でもっと協力ができたかもしれないので勉強したいとこですね。またこういう機会があればぜひ参加したいです。

2016年 1月
単位を取るために直前に必死こいて勉強したりしてた。
特筆するべきことはなし。

2016年 2月〜3月
バイトをして金を稼いだ。
朝永ゼミからサクライゼミにいつの間にかシフトしていた。(ブラケット記法は素晴らしい)
時折遊んだ。


こうして去年一年を振り返ってみるとなにかしたと言えるのはHackUにでたことぐらいでしたね。

今年はもう少しいろいろ手を出したいところですね。
まずは研究室選ばないと。

*1:今年も計算機コースの歓迎会は開かれない模様です。にんまり。

Stokesの定理(一般)

今日勉強した非常に美しい定理

\(\omega\)を多様体\(\mathcal{M}\)上の\(p\) 形式、\(c\)を\(\mathcal{M}\)上の\( p+1 \) 鎖とすると次が成り立つ。


\[
\int_{\partial c} \omega=\int_c d\omega
\]


簡潔すぎてかえってわかりにくいw

極座標ラプラシアンを求める公式

座標として \( {\bf x}(q_1,q_2,q_3)\) が入っているとき

\[
\newcommand{\partdif}[2]{\frac{\partial #1}{\partial#2}}
\]

\[\partdif{{\bf x}}{q_1} ,\partdif{{\bf x}}{q_2}, \partdif{{\bf x}}{q_3} \]が

各点で互いに直交しているとき、ラプラシアンを求めるのに次の公式が使えます。




\[
\lambda_1{\bf e}_1=\partdif{{\bf x}}{q_1},\lambda_2{\bf e}_2=\partdif{{\bf x}}{q_2},\lambda_3{\bf e}_3=\partdif{{\bf x}}{q_3}
\]

であるとし

\({\bf e}_1,{\bf e}_2,{\bf e}_3\)が右手系正規直交基底となるように取ると

\[
\Delta f=
\frac{1}{\lambda_1\lambda_2\lambda_3}\left[\partdif{}{q_1}\left(\frac{\lambda_2\lambda_3}{\lambda_1}\partdif{f}{q_1}\right)+\partdif{}{q_2}\left(\frac{\lambda_3\lambda_1}{\lambda_2}\partdif{f}{q_2}\right)+\partdif{}{q_3} \left(\frac{\lambda_1\lambda_2}{\lambda_3}\partdif{f}{q_3}\right)\right]
\]






導出はひとまず置いておいて、実際に使ってみます。

極座標

\[
\begin{cases}
x_1=r\sin \phi \cos \theta\\
x_2=r\sin \phi \sin \theta\\
x_3=r\cos \phi \\
\end{cases}
\]

で導入すると
\[
\begin{cases}
\partdif{{\bf x}}{r}=(\sin \phi \cos \theta,\sin \phi \sin \theta,\cos \phi ) \\
\partdif{{\bf x}}{\phi}=(r\cos \phi \cos \theta,r\cos \phi \sin \theta,-r\sin \phi) \\
\partdif{{\bf x}}{\theta}=(-r\cos \phi \sin \theta,r\cos \phi \cos \theta,0) \\
\end{cases}
\]

\[
\begin{cases}
{\bf e}_1=(\sin \phi \cos \theta,\sin \phi \sin \theta,\cos \phi ) \\
{\bf e}_2=(\cos \phi \cos \theta,\cos \phi \sin \theta,-\sin \phi) \\
{\bf e}_3=(-\sin \theta,\cos \theta,0) \\
\end{cases}
\]

と置くと

\[
\begin{cases}
\lambda_1=1\\
\lambda_2=r\\
\lambda_3=r\sin \phi\\
\end{cases}
\]

と置けばよく、公式に代入すると


\[
\Delta f=
\frac{1}{r^2\sin\phi}\left[\partdif{}{r}\left(r^2 \sin\phi\partdif{f}{r}\right)+\partdif{}{\phi}\left(\sin\phi\partdif{f}{\phi}\right)+\partdif{}{\theta} \left(\frac{1}{\sin\theta}\partdif{f}{\theta}\right)\right]
\]

線形微分方程式の解の次元について

とりあえず、微分方程式の解の存在と一意性の定理を使えば
\(n\)次線形微分方程式の解空間の次元が\(n\)であることはすぐ示せるとわかった。

こんなかんじでやる。

線形微分方程式

\[
\frac{d^{n}y}{dt^{n}}+a_1\frac{d^{n}y}{dt^{n}}+\cdots+a_ny=0
\]

が初期条件

\[
\frac{d^{j}y}{dt^{j}} =b_j(t=t_0)
\]

\[
(j=0,\cdots,n-1)
\]


を満たす場合を考える。

解の存在性と一意性からこの条件を満たす解\(y^\circ \)がただ一つある。

いま解の存在性から

\[
\frac{d^{j}y_i}{dt^{j}}=
\begin{cases}
1&&(i=j)\\
0&&(i\neq j)
\end{cases}
(i,j=0,\cdots,n-1)
\]

となるような線形独立な解の組\(y_1,\cdots,y_n\)

がとれて

\[
y^\dagger=b_1y_1+\cdots+b_ny_n
\]

と置けばこれは

微分方程式と初期条件を満たす解になっている

よって一意性より

\[
y^\circ=y^\dagger
\]

となり、\(y_1,\cdots,y_n\)は解空間を張る。

よって解空間の次元は\(n\)となる。

').size(unit * r, unit * r).after(loop); loopFlag = !loopFlag; } focusCircle.stop().animate(450, '>', 250).opacity(0.5).after(loop); pressedCircle.stop().animate(200).size(unit * radius, unit * radius).opacity(0); } function onBlur() { focusCircle.stop().animate(800).opacity(0).size(unit * radius, unit * radius); } function onPress(event) { var innerRadius = radius - 12; var outerRadius = radius + 4; pressedCircle.center(event.layerX, event.layerY).opacity(0.2).size(unit * innerRadius, unit * innerRadius) .animate(600, '>').size(unit * outerRadius, unit * outerRadius).opacity(0); focusCircle.stop().animate(200).opacity(0).size(unit * radius, unit * radius); } // イベント設定 group.click(onPress); group.touchstart(onFocus); group.touchend(onBlur); group.mouseover(onFocus); group.mouseout(onBlur); />