スパースモデリングとは、スパース性を用いた統計的データ分析手法である。様々な機械学習の手法と組み合わせることで驚異的な威力を発揮する。
概要
スパースモデリングは機械学習法の一種で、「大量のデータの中から小数の本質的なデータを高い確率で抽出する方法」である。まばらであることをスパース(sparse)といい、本当に重要な情報は全体のほんの僅かでまばらにしか存在しないという仮定に基づいている。
AIでおなじみのディープラーニングとの違いは、ディープラーニングは非常に「それっぽい」データを吐きだすがその機序はよくわからない一方、スパースモデリングはブラックボックス的なアルゴリズムが無く、何をやっているか見通しが立てやすいがスパース性という仮定に基づいている、という点であろう。うまく「少数の本質的な情報」を抜き出せるようにデータの前処理や事前の辞書学習、アルゴリズムの工夫が必要である。
ディープラーニング同様アルゴリズムの工夫次第で、ノイズ除去、データ圧縮、特徴量抽出、欠損部分の推定、学習無し処理など、様々なことを柔軟にこなすことができる。
劣決定逆問題
線形代数の最も基本的な方程式は以下の式で表される。
y=ax
aはある定数であり、データxが与えられた時のyの値を求めよ、という式である。xが実数の時、yも実数であり、解はただ一つに定まる。
逆に、yが与えられた時、xは何であったか?ということを調べたい時がある。これを逆問題と呼ぶ。答えはx=y/aと変形することで求まる。
現実では、「牛丼一杯a[kcal]である。今、牛丼をx[杯]食べた。何[kcal]摂取したか?」→逆問題化→「牛丼一杯a[kcal]である。今、y[kcal]摂取した。牛丼は何[杯]食べたか?」という問題などが考えられる。
xとyが1次元の時は以上おしまい、となるが、xとyがベクトルの時は次元次第で解がただ一つに決まったり決まらなかったりする。
xとyの次元が等しい時、aは正方行列で表現され、行列aが正則であれば逆問題を解くことができる。しかし、xの次元>yの次元であるとy=axの逆問題を一意に解くことは不可能となる。例えば、2個のデータ(y)から3個のデータ(x)を一意に復元することはできない。3次元空間の点xを表現するのに3個の基底ベクトルが必要なのに、yからは最大2個のベクトルしか生成できず、足りない分不確定になるためである。
このように、得られた点yのデータ数がxを復元するのに不十分な数であるとき、劣決定逆問題であるという。
現実では画像や音源の劣化、データの不可逆圧縮などで見られる現象である。
スパース性
スパースとは、粗であること、スカスカであることを言う。例えば、10000点の0と1からなるデータxの内、100個が1で残りは全て0であるとき、データxはスパースであるという。
音声や画像を周波数空間で表現するとスパースになりやすいことが知られている。
スパースな劣決定逆問題
先程の劣決定性は、xの次元がyの次元より大きいために起こる。しかし、xがスパース性を持つ時、解が一意に求まることがある。例えばxが3次元の縦ベクトルの時、以下の逆問題
Y=AX
|
= |
|
|
は2個の方程式で3個の未知数を求めるという問題になり、これでは解が一意に決まらない。しかし、「xの非ゼロ要素が一つだけである」と分かっているならば、(x1 x2 x3)=(0 0 1)と一意に決定される。
このように、元となるデータXがスパースであるという仮定の下では劣決定性があっても逆問題を解くことが可能になる場合がある。
スパースモデリングとは、「スパース性を仮定することで、得られた観測データから元のデータを復元する手法」と言うことができる。しかし、非ゼロ要素をあてずっぽうに探し出すのは非常に計算コストがかかる困難な問題であり、そのような問題を解く効率の良い近似アルゴリズムがスパースモデリングということである。
事後確率推定
劣決定逆問題は、「観測Aの下で観測データYが得られたとき、元のデータXはどのようなものであったか?」という問題を解くことである。推定されたデータXcは元のデータXと一致する確率が他の候補と比べると最も高いことが期待される。これは事後確率最大推定と呼ばれる。スパースモデリングのアルゴリズムはこの事後確率を最大化するものが中心である。
Y=AXの式で表されるYを生成するデータXcは無数にある。また、観測には何かしらの誤差εがつきものである。その中で最もよく再現するものは、YとAXcの差(Y-AXc)=εが十分小さいことが期待される。そこで、ある非負実数値を返す「ペナルティ関数」fを考え、f(Y-AXc)つまり誤差が最も小さくなるようなXcを選ぶという戦略が考えられる。ペナルティ関数には通常、各成分の二乗和(つまり普通のユークリッド距離関数の2乗)が使われる。
しかし、Xcが非現実的なものであっては困るので、f(Y-AXc)がただ小さければよいというものではない。そこで、新しいパラメータ関数Sを用意する。Sは必ず満たしてほしい条件(最小、最大、特定の値、など)を記述する関数である。すると元の問題は以下のように書き直される。
min[ f(Y-AXc) ] ただし S(X)が特定の条件を満たす
これは条件Sを満たしつつペナルティの最も小さいXcが求めたいXだ、という式である。
ラグランジュの未定定数法を使うと、実数λを用いることで一つの式で書けるようになる(S(X)が小さいほど良い近似となるように工夫しておく)。
この式は「Xcが得られた信号にどれだけ近いか」という成分と「Xcが要求される尤もらしさをどれだけ持っているか」という成分をうまく混合しているという式である。λが小さければ前者を優先したデータとなり、λが大きければ後者を優先したデータとなる。
関数Sがエントロピーと呼ばれるものである場合、最大エントロピー法になる。言い換えれば、最大エントロピー法はエントロピーが最大という条件の下で解く劣決定な事後確率最大推定法ということになる。最大エントロピー法の弱点は、非ゼロ成分がノイズのように大量に含まれる傾向にあり重要な量だが絶対値は小さいような成分が隠れてしまうため、必ずしも「エントロピー最大→最もよい近似データ」とはならないことにあった。
l0ノルム
ベクトルを引数に持ち非負の実数値を返す関数をノルムと呼び、そのようなノルムを備えたベクトル空間のことをノルム空間と呼ぶ。典型的なノルムはユークリッドノルムl2(X)=|X|2=∑i(xi2)である。この指数2は非負実数に拡張することができ、実数pであるときのノルム∑i(xip)をlpノルムという。
また、ベクトルの非ゼロの要素数を返すノルムをl0ノルムという(∑i(xi0)という形はしていないがp→0の極限でlpノルムはl0ノルムになる。)
LASSO回帰
ペナルティ関数Sがノルムlpであるとき、解きたい式は以下になる。
スパースモデリングは「現実を表すデータ数は少ない」つまり非ゼロ要素が非常に少ない(スパースである)という仮定に基づく。lpノルムは、pが1以下であるとき、Xがスパースとなることが保証されている。スパース性が保証される。lpノルムによりスパースなXcを求めるアルゴリズムを圧縮センシングという。
関数Sとしては非ゼロ要素をそのまま返すl0ノルムが最良である。しかし、どれが非ゼロかは全ての組み合わせを総当たりで調べるしかない。これは多項式時間内に終わらない非常に時間のかかる問題となる。なので、l0ノルムの代わりにl1ノルムを使った方式が採用され、これをLASSO(ラッソ、Least Absolute Shrinkage and Selection Operators)という。
実装
詳しい解説は専門書に任せるとして、LASSO回帰のアルゴリズムの一つであるADMM for LASSOのScilabによる実装例を記述する。これはコードが短いわりに収束が早く、しかも不安定性がないといわれる優秀なアルゴリズムである。
まず、スパースな原信号X0をそのままXの形で観測し、観測の際に加法的ノイズεが乗るという以下のモデルを考える。
Y=X0+εであり、原信号を復元するデータXを反復回数tにより求める。mu=0.3、lm=0.3、t=100程度でノイズが大まかに除去できることがわかる。
人工的に測定ノイズを含むデータを作成しただけなので、もしYとして実測したデータを使えるならノイズεを足す必要はない。
//=======================================メイン=========================
function [x,y,z,x0,u]=admm1(mu,lm,t)
n=1000 //データサイズ
k=10 //非ゼロデータ数
x0=[rand(1,n,"normal") zeros(1,n-k)]
x0([1:1:n])=x0(grand(1,"prm",(1:n)))
y=x0'+rand(n,1,"normal")/10 //非ゼロ成分+加法的白色ノイズ
x=zeros(n,1)
z=zeros(n,1)
u=zeros(n,1)
for i=1:t
x=y/(1+mu*lm)+(z-u)*(mu*lm/(1+mu*lm))
z=SoftThr(x+u,1/mu)
u=u+x-z
end
endfunction
//=======================================軟判定閾値関数=========================
function y=SoftThr(x,l)
y=x
y=zeros(y)
temp=find(x>l)
if temp<>[] then
y(temp)=x(temp)-l
end
temp=find(x<-l)
if temp<>[] then
y(temp)=x(temp)+l
end
endfunction
次に、周波数空間でスパースな2n点のランダムデータX0を観測し、n点周期データYを測定した場合を考える。
Aはフーリエ変換を表す観測行列であり、Y=AX0+εである。復元したい周波数信号Xが変換Aを介してYと十分に近く、さらにXはスパースに表現できるという設定である。こちらもmu=0.3、lm=0.3、t=100程度でノイズを大まかに除去できる。
下記の観測行列Aは、フーリエ変換行列をcos成分とsin成分に分けてn×2n行列として表現している。観測行列Aがn×nの単位行列である場合、先ほどの例と一致する。
観測行列はフーリエ行列を始め、離散コサイン、ウェーブレットなどの各種基底セットや、サンプルデータから学習した辞書などが用いられる。
//=======================================メイン=========================
function [x,y,z,u]=admm(mu,lm,t)
n=1000
k=10
nn=[0.000001:1:n]
A=[cos(2*%pi*nn'*nn/n) sin(2*%pi*nn'*nn/n)]
x0=[rand(1,n,"normal") zeros(1,n*2-k)]
x0([1:1:n*2])=x0(grand(1,"prm",(1:2*n)))
y=A*x0'+rand(n,1,"normal")/10
A1=inv(A'*A/lm+eye(2*n,2*n)*mu)
A2=A1*A'/lm
x=zeros(2*n,1)
z=zeros(2*n,1)
u=zeros(2*n,1)
for i=1:t
x=A2*y/lm+A1*(z-u)*mu
z=SoftThr(x+u,1/mu)
u=u+x-z
end
endfunction
関連項目
- 2
- 0pt