読者です 読者をやめる 読者になる 読者になる

北の大地から

北大でひっそりと研究

NMFアルゴリズムの紹介　その② GCD

研究が進んでいないので，記事を書いて，モチベーションを高めることに．．．．
今回の紹介アルゴリズムはNMFで現在(2015年10月)で最速といわれている(筆者の知る限り）手法Greedy Coordinate Descent Algorithmの紹介です．
元論文は，
Hsieh, C.-J., & Dhillon, I. S. (2011). Fast coordinate descent methods with variable selection for non-negative matrix factorization. In ACM SIGKDD (pp. 1064–1072). https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=3&cad=rja&uact=8&ved=0CDUQFjACahUKEwix5Oar-arIAhVJmZQKHcJJBik&url=https%3A%2F%2Fwww.cs.purdue.edu%2Fhomes%2Fdgleich%2Fconf%2Fslides%2Fmmm2012%2Fdhillon.pdf&usg=AFQjCNGssbMCO2vHvbaOIl3nkDCnuE6JiA&sig2=Tti3GHlrtBTuYGf8bVbAhA&bvm=bv.104317490,d.dGo
Dhillon先生は，超やばい研究者で，そのお弟子さんたちと合わせて，毎年恐ろしい本数の論文をトップ会議へ通してます．

それはおいといて， Greedy Coordinate Descentは要素毎に更新する手法です．
HALSも要素毎といえば要素毎の更新といえばそうなのですが，
何故，敢えて区別していうかは，読んでいただければわかってくると思います．
GCDの肝は，更新する要素を選択して重要な要素をgreedyに更新を行うところです．
HALSもMUも一回の更新で全ての要素を平等に更新します．
一方，GCDでは近似精度を上げるために重要な要素はそれが重要である限り何度でも更新し，不要であれば一切は更新しません．
そうすることで劇的に収束が早くなっています．

* GCD
では実際に更新方法を導出してきます．
その前に，どのような方針で最適化するかを大まかに述べます．（更新する行列はUとします）

  1. 行列Uの全ての要素に対して，要素毎の問題の最適解を計算
  2. 最適解を用いて各要素の重要度を計算 +一番重要な要素を更新
  3. 更新された要素に関係のある要素の最適解及び重要度を更新
  4. 2-3を繰り返す

これだけです．重要なのはStep3が実は高速に行えるという点です．

導出

まず，毎度お馴染みのNMF問題は
\begin{equation} \| X - UV^{T}\|^{2}_{F} \end{equation}
です．ちなみに転置するとUもVも同じなので，今回は本当にUのみの話をします．

要素毎の最適解を見つけるのでU_{ij}についての最適化を考えます．（ X_{ij}毎の問題ではありませんので注意．）
最適なU_{ij}の値（他を固定している場合）はHALSのものと同じ結果になります．（ここが下線を引いた部分に対応します．どうしてそうなるのかを考えると，GCDが高速に分解できる理由が見えてきます．）

具体的には，GCDの論文の書き方にならって書くと，
\begin{equation} U_{ij} \leftarrow max(0, U_{ij} - (\frac{(UV^{T}V-XV)_{ij}}{(V^{T}V)_{jj}})). \end{equation}

さらに，更新の差分であるsは，
\begin{equation} s=U_{ij}^* - U_{ij} = max(0, U_{ij} - (\frac{(UV^{T}V-XV)_{ij}}{(V^{T}V)_{jj}})) -U_{ij}. \end{equation}
と書けます．

あっさりではありますが，Step1は突破しました．

次にStep2です．GCD手法の最大の貢献はこの部分をテイラー展開であっさりと計算できるようにしたところです．（他分野ではあったと思いますが，NMFにはここで初めて持ち込まれました．）
U_{ij}を上の更新式で更新した場合の近似誤差の関数はsを引数とした関数g_{ij}(s)で書けます．
\begin{equation} g_{ij}(s)=\frac{1}{2}\sum_{k}(X_{ik}-(UV)_{ik}-sV_{kj})^{2} \end{equation}

この式を2次以降は無視したテイラー展開で近似します．
\begin{equation} g_{ij}(s)=g_{ij}(0)+g_{ij}'(0)s+\frac{1}{2}g_{ij}''(0)s^{2} \end{equation}
そうすると，どれだけ目的関数を減らすかということが，差分で分かるようになります．その重要度 D_{ij}とすると
\begin{equation} D_{ij}=g_{ij}(0)-g_{ij}(s)=-g_{ij}'(0)s-\frac{1}{2}g_{ij}''(0)s^{2} \end{equation}

ここで，g_{ij}'(0)とg_{ij}''(0)の値はどうなっているかというと，

\begin{equation} g_{ij}'(0)=(UV^{T}V-XV)_{ij} \end{equation}
\begin{equation} g_{ij}''(0)=(V^{T}V)_{jj} \end{equation}
です．


ここで，GCDが高速に計算できる理由は，行列Uの更新において，互いに影響する要素は同じ行ベクトルに存在する要素のみということです．
具体的には， U_{ij}を更新したと時に，再計算が必要なのは\forall \ k \ \ \ g_{ik}'だけということです．

実際の更新は
\begin{equation} \forall \ k \ \ g_{ik}'(0)=g_{ij}'+s(V^{T}V)_{jk} \end{equation}
です．

行ベクトルのみしか影響しないため，Step2とStep3は行ベクトル内で行います．
まとめると

行列Uの更新は

  1. A=UV^{T}V-XV
  2. B=V^{T}V
  3. S_{ij}=\max(U_{ij}-\frac{A_{ij}}{B_{jj}},0)-U_{ij}
  4. D_{ij}=-A_{ij}S_{ij}-\frac{1}{2}B_{jj}S_{ij}^{2}
  5. 各行ベクトルに対して(i行目を載せる)
  6. while convergence
  7. k=arg\max_{j}D_{ij}
  8. U_{ik} \leftarrow U_{ik}+S_{ik}
  9. A_{i:} \leftarrow A_{i:}+S_{ik}B_{k:}
  10. S_{ij}=\max(U_{ij}-\frac{A_{ij}}{B_{jj}},0)-U_{ij} 全てのjに対して．
  11. D_{ij}=-A_{ij}S_{ij}-\frac{1}{2}B_{jj}S_{ij}^{2}　 全てのjに対して．
  12. while end

行ベクトルでのループの更新は全てO(k)で高速にできるため，高速に計算できます．

実験

実験の比較です．
前々回の記事で，紹介していますが，GCDのコードは著者によって提供されています．www.cs.utexas.edu

注意して欲しいのは，GCDはfor文を多用することになるので，MATLABやRでの実装は向きません．
実際に，提供されているコードはMATLABでCを呼び出すmexを用いています．（厳密に言えば公正な比較ではありません）
それでも，尚計算量の解析は元論文でされているのでそちらを参照してください．

ちなみに更新に提供されている関数は doiterと呼ばれるものです．
使い方は　
mex doiter.cpp
しといて，
U^{T} = doiter(A^{T},B,U^{T},tol,K^{2})
という形で用います．後ろ二つは更新の打ち切りに用いるものなのであまり影響を受けません．おまじないと思ってください．
詳しい話は論文に載っています．
tol=0.001として使って大丈夫です．

比較手法

MUとHALSの比較を行います．
MUについては，他のサイトの記事，HALSについては前回の記事を読んでください．

評価指標

前回と同じです．

  • CPU時間
  • 近似精度 Normalized Residual Value

データセット

前回と同様にReuters21578データセットです．
X \in \mathbb{R}^{18993 \times 8293}の行列を分解します．

実験設定

同様に分解後のサイズKを変動させて挙動を観測します．
いくつか Kを固定した結果も載せます．

実験結果

まずは，全てのK=10,\dots100において，100回の反復更新を行った場合の結果です．
f:id:Kei210:20151005195335p:plain
HALSがK=100のときにのみ急激に上がるのが気になりますね．．．
まぁおいといて，一回の計算量は何度も要素を更新している点で，遅いことが分かります．

K=30の際の実際の計算時間はこちらです．
f:id:Kei210:20151005195512p:plain

さて次は，更新によってどれだけ近似誤差が減少しているかです．
K=50のときは
f:id:Kei210:20151005195908p:plain

・・・あれ？
GCDがそんなに良くない．．．おかしいな．．．．

計算時間でも
f:id:Kei210:20151005195722p:plain

HALSに負けてる．．．？
ちょっと原因を探している状態です．．．　おかしいな．．．
分かった方は教えてください．．．

まとめ

GCDは要素毎の更新を行うNMFです．
各要素に対して，重要度を算出して，それに沿ってgreedyに更新を行う手法です．
同一の行列内では，同じ行indexを持つ要素のみしか，互いに影響しないため，行ベクトル毎にループをするだけでよく，
これにより，高速に要素の更新＋重要度の更新を行うことができます．

実験では高速化されるはずだったのですがされてませんでした．失敗している模様です．

コードは以下からダウンロードできます．前回のリンクと同じですがアップデートした形になります．
ただdoiter関数は含めていないので，使う場合は，著者のページからダウンロードして用いてください．
https://dl.dropboxusercontent.com/u/97469461/NMF.zip


NMFGCDのコードです．

function[U,V]=NMFGCD(X,U,V,iter);

%%Input
%% X = M x N matrix
%% U = M x K matrix
%% V = N x K matrix
%% iter: #iterations


%% initialization 
K=size(U,2);
eps=1e-08;



for loop=1:iter

	
	B=V'*V;
	A=(U *B) - (X*V);
	Unew=doiter(A',B,U',0.001,K^2);
	U=Unew';
	
	B=U'*U;
	A=(V*B) -(X'*U);
	Vnew=doiter(A',B,V',0.001,K^2);
	V=Vnew';
	
end