# [論文紹介]変数選択を用いた高速な非負値行列因子分解 + α

[ブログ](https://localmin.github.io/2018/12/11/fast-nmf/)の方で数式を扱えるようになったので加筆修正は今後は[ブログ](https://localmin.github.io/2018/12/11/fast-nmf/)で行う予定です。こっちの更新は行いません。悪しからず。
まあ今のところ大きな変更はないですが。

## はじめに
おばんです！ [localmin](https://twitter.com/localmin?lang=ja)と申します。現在M2で、普段の研究ではあえてぼかした言い方をすると行列演算について考えたり、パソコンをカタカタしています。本記事は[#今年読んだ一番好きな論文2018](https://adventar.org/calendars/3392)の11日目の記事です。本記事でも少し古いですが行列にまつわる論文について紹介していきたいと思います。

## 紹介論文[1]
* 題名
[Fast Coordinate Descent Methods with Variable Selection for Non-negative Matrix Factorization](http://www.cs.utexas.edu/~cjhsieh/nmf_kdd11.pdf)
* 著者
Cho-Jui Hsieh, Inderjit S. Dhillon
* 出典
Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining Pages 1064-1072, 2011

## 論文概要
非負値行列因子分解(Non-negative Matrix Factorization)(以下NMFとします)の高速化に関する論文です。NMFは反復更新式を用いて最適化を行いますが、既存手法では全変数を等しく更新(要するに更新しなくてもよい要素を更新してしまう)という欠点があります。この論文では目的関数にフロベニウスノルムの二乗を用いたときに、必要な変数のみを選択し更新を行う手法を提案をしています。
また目的関数にKL-divergenceを用いた場合の新しい手法も提案してます(ただし変数選択は用いない)。
各提案手法と既存手法を実データを用いて比較を行い、提案手法の収束の速さを確認してます。

...と小難しく書きましたが、これからなるべくわかりやすく説明していく予定です。たぶん。

## なぜこの論文を選んだか？
もちろん今年読んだ中で印象に残った論文だからというのもありますが、
1. 巷に溢れるNMFの記事には基本的にMultiplicative Updateの記事しかないから(然るべき理由はありそう)
2. 速いアルゴリズムはかっこいいから

といった理由もあります。この記事をきっかけにNMFにはMultiplicative Updateだけでなく、速いアルゴリズムがあるよ！ということを少しでも理解してもらえば冥利に尽きます。
っていうか速いアルゴリズムってかっこいいじゃないですか！！(個人の感想)

## 非負値行列因子分解(NMF)[2]とは
NMFについては論文内でも説明されていますが、少しあっさり風味なのでここでは少し詳しめに説明します。
世の中のデータには非負値の行列で表現できるものが多々あり(例えばテキスト、画像、音響信号のパワースペクトルなど)、その特徴を解析したいといった需要があります。
NMFは各要素が非負で与えられる行列を二つの非負値行列の積に近似分解してあげる手法です。具体的には各要素が非負である観測行列$V\in {\bf R}^{m \times n}$を以下のように近似分解します。

$$
V \approx WH
$$

$W\in {\bf R}^{m \times k}$は基底行列、$H\in {\bf R}^{k \times n}$は係数行列となっています。
NMFはいわゆる次元削減と呼ばれる手法に分類され、観測行列をより小さなランクの行列で近似することを目的とします。
$k$は近似に用いる基底ベクトルの本数をあらわしており、何本の基底ベクトルで近似するかは使用者が事前に指定します。
行列に馴染みがある人は基底行列、係数行列と言われてすぐにわかるかもしれませんが、そうでない人のためにもう少し詳しく説明すると、
観測行列$V$の$i$列目の列ベクトルを$v_{i}$、$H$のj列目の列ベクトル$h_{j}$とすると、NMFでは観測行列の各列は

$$
v_{i} \approx Wh_{j}
$$

と表現することができます。
すなわちNMFは観測行列の各列を$W$の列ベクトルの線形結合で近似していると解釈することができます。
線形結合に使用されるベクトル(基底)を格納した行列なので基底行列、そして線形結合に置ける係数を格納した列ベクトル$h_{j}$を並べているので係数行列というわけです。

NMF以外にも行列を分解する手法で有名なものに特異値分解(SVD)や固有分解(スペクトル分解)などがあります。しかし非負値行列を適用した場合に、分解した結果に負の値が出現してしまったり、得られるベクトルが直交してしまうため、結果の解釈が難しくなる欠点が存在します。
そこでNMFでは与えられたデータはもちろん、分解結果にも非負の制約与えることでデータの見通しが良くなる長所があります。また分解して得られた基底ベクトルは必ずしも直交しません。

## NMFの解法
### 目的関数の設定
次にどうやって$W$、$H$を求めるのかという話に移っていきます。
人によっては行列を二つの積にバラすのは簡単では？と思う人もいるかもしれません。でもそれは具体例をみれば難しいとわかるでしょう。
では以下の方程式をを$k=2$という条件で二つの非負値行列$W$、$H$を求めてみてください。

$$
\begin{pmatrix}
23 & 13 & 22\\
8 & 16 &40\\
8 & 5 & 9
\end{pmatrix}=WH
$$

どうですか？難しくないですか？っていうかそもそも解あるのか？ってなると思います。実際のデータは何百、何千のサイズになるのでもっと大変です。

実際にこのような行列 $W, H$を求めることは可能とは限りません。しかし不可能であっても近似的に求めることは可能なはずです。

ここで目的関数を用意します
* フロべニウスノルムの二乗

$$
f(W, H) = \frac{1}{2}||V - WH||_{F}^{2}
$$

$||\cdot||_{F}$ はフロべニウスノルムと呼ばれ、行列の全要素の二乗和の平方根をとります。

この式は観測行列$V$と分解して得られる行列$WH$の距離を測る関数になっています。すなわち分解して得られた$WH$が観測行列に近いほど上述の関数の値が0に近づきます。
よってNMFは上述の関数を最小化するような$W,H$を求める問題になります。残念なことにそのような$V, H$は解析的に求めることができないため反復更新によって近似的に求めます。

### 一変数ごとの座標降下法
この論文では反復更新の方法として一変数ごとの座標降下法を用いています。具体的には

$$
(W,H)\leftarrow (W + sE_{i,r}, H)
$$

のような更新を収束するまで行います。$E_{i,r}$は$(i,r)$成分が1, 他の成分が0である$m\times k$の行列です。
座標降下法ではある一つの要素に注目して他の変数は固定し、その変数での部分最適化問題を解くということを繰り返していくことになります。
この式から読み取れるのは$W$の$(i,r)$成分(以下$W_{i,r}$)を$s$だけ更新しているということ、すなわちまず僕らが知りたいのは$s$の具体的な値になります。

ではどんな$s$が欲しいのか、それは目的関数を減少する$s$に他なりません。そのような$s$を求めるためには以下のような一変数の部分問題を解きます。

$$
\min_{s:W_{i,r}+s\geq 0}g_{i,r}^{W}(s)\equiv f(W+sE_{i,r}, H)
$$

要するに$W_{i,r}$を$s$だけ更新した時、目的関数を最小化するような$s$を求めることになります。
$g_{i,r}^{W}(s)$を$s$について整理してあげると、

$$
g_{i,r}^{W}(s) = g_{i,r}^{W}(0) + (g_{i,r}^{W})^{’}(0)s + \frac{1}{2}(g_{i,r}^{W})^{’’}(0)s^2
$$

となり、$s$に関する二次関数になります。こうなってしまえば簡単で、$s$について微分し、その導関数を0にするような$s$を求めます。ただし、$W_{i,r}+s\geq 0$なる制約があるので

$$
s = \max \left(0, -\frac{(WHH^{T}-VH)_{i,r}}{(HH^{T})_{r,r}}\right) - W_{i,r}
$$

という結果になります。max関数を使用することで非負性を保持できます(実際に値が負になるような数値を代入してもmax関数のおかげで0になります)。
これによって

$$
W_{i,r}\leftarrow W_{i,r} + s
$$

と更新することができます(ちなみに $V$についても同様)。

## 変数選択による高速化
### Outer Iteration & Inner Update
従来手法(Multiplicative Update[2], FastHALS[4])では上に似たような更新式を導出した後、全変数に対して等しい回数更新を行います。しかし、その場合更新しなくても良い変数を更新してしまうという欠点があります。ここでいう更新しなくても良い変数とは、更新しても目的関数の値が減少しない変数のことであり、そのような変数を更新してしまうことで余計な計算量をくうことになります。
よってこの論文では重要な変数、すなわち目的関数の値を大きく減少させることができる変数のみを選択し貪欲に更新し続ける手法を提案しています。

まず変数選択法の説明の前に更新手順について説明します。この論文では便宜上、注目する行列を変えることをOuter Iteration、注目している行列内での更新をInner Updateとしています。
Outer Iterationでは$W, H$を交互に更新していきます。

$$
(W^{0}, H^{0}) \rightarrow (W^{1}, H^{0}) \rightarrow (W^{1}, H^{1})...
$$

そして一回のOuter Iterationの間にInner Updateという行列内の更新を行います。

$$
(W^{t,0}, H^{t-1}) \rightarrow (W^{t,1}, H^{t-1}) \rightarrow (W^{t,2}, H^{t-1})...
$$

実際の変数選択はInner Update内で行います。

要するに、
1. $W$の各要素を変数選択を用いて十分に更新を行う($W$のInner Update)
2. $H$の各要素を変数選択を用いて十分に更新を行う($H$のInner Update)
3. また$W$に戻って同じことを繰り返す

といった処理を目的関数の値が収束するまで繰り返す、これが基本的な流れになります。

### 変数選択の方法
$W$の変数選択を行う場合を考えます($V$についても同様に求められます)。
変数選択の方法はシンプルで、目的関数の減少量が多い変数を選択します。そのために$W$について、各要素を更新した際の目的関数の減少量を格納した行列を$D^{W}$を用意します。サイズは$W$と同じで、この$(i,r)$成分には$W_{i,r}$を更新した時の目的関数の減少量が格納されています。
$D^{W}$の各成分(目的関数の減少量)以下のように計算されます。

$$
D_{i,r}^{W} \equiv g_{i,r}^{W}(0) - g_{i,r}^{W}(s)
= -G_{i,r}^{W}s-\frac{1}{2}(HH^{T})_{rr}s^{2}
$$

$g_{i,r}^{W}(s)$は$W_{i,r}$を$s$だけ変化させたときの目的関数の値なので上式が更新後の目的関数の減少量を計算していることはわかると思います。(ここでは$G^{W}=WHH^{T}-VH^{T}$)。
当然のことながら更新を行えば$D^{W}$の値も変化していきます。なので$D^{W}$を常に最新に保つことができれば、これに基づいて変数選択を行うことができます。しかし更新の度に先ほどの式に従って$D^{W}$の値を全て再計算していると$\mathcal{O}(mk)$ かかってしまいます。たとえどんなに更新一回あたりの目的関数の減少効率がよくても、変数選択にに時間がかかっていては本末転倒です。
よってこのオーバーヘッドを軽減するために$W_{i,r}$を更新した時、$G^{W}$の$i$行目しか変化しないという特徴を利用します。この特徴から$D^{W}$も$i$行目しか変化せず、再計算する必要があるのは$G^{W}$の$i$行目だけになります。
以下の式で$G_{i,j}^{W}$を更新することで$D^{W}$を最新に保つことができます。

$$
G_{i,j}^{W}\leftarrow G_{i,j}^{W} + s(HH^{T})_{r,j} \forall j=1,\dots,k
$$

これらの計算量は$\mathcal{O}(k)$であり、$D^{W}$全ての要素を再計算するよりも少ないコストで変数選択をすることができます。
これがこの論文のキモであり、変数選択による計算の軽減に成功しています。

このことから更新が各行列に及ぼす影響は$i$行目だということがわかります。なので基本的には、行ごとに変数選択および更新を行えば良いということになります。
そのため$i$行目の更新において次に選択される列インデックス$q_{i}$は

$$
q_{i} = arg\max_{r}D^{W}_{i,r}
$$

のようにして選択されることになります。

このように変数選択しながらを更新していく手法をGreedy Coordinate Descent(以下GCD)と論文では名付けています。いいかんじですね。

## 実験
いろいろややこしいことを書いてきましたが要は速いのかって話ですよ。実データを用いて実験をします。この論文では二種類の実験を行なっています。
一つは密なデータに対して、もう一つはスパースなデータに対してです。

### 実験1(密なデータに対して)
比較する手法はGCD、FastHALS[4]、ProjGrad[5], BlockPivot[6]（あれ？Multiplicativeは？)。
各アルゴリズムにおいてデータセット毎に設定した相対誤差になるまでの計算時間及び浮動小数点数演算の回数(FLOPs)を比較します。
もちいたデータは以下のようになっています。
1. Synth03 : 乱数によるデータ。30％が0を要素の持つ
2. Synth08 : 乱数によるデータ。80％が0を要素の持つ
3. CBCL : 画像データ
4. ORL : 画像データ

|dataset|m |n | k |relative error|GCD|FastHALS|ProjGrad|BlockPivot |
|---|---|---|---|---|--- |---|---|---|
|Synth03|500|1000 | 10 |$10^{-4}$ |0.6/0.7G | 2.3/2.9G | 2.1/1.4G | 1.7/1.1G |
|Synth03|500|1000 | 30 |$10^{-4}$ |4.0/5.0G | 9.3/16.1G | 26.6/23.5G | 12.4/8.7G |
|Synth08|500|1000 | 10 |$10^{-4}$ |0.21/0.11G | 0.43/0.38G | 0.53/0.41G | 0.56/0.35G |
|Synth08|500|1000 | 30 |$10^{-4}$ |0.43/0.46G | 0.77/1.71G | 2.54/2.70G | 2.86/1.43G |
|CBCL|361 |2429 | 49 |0.0410 |2.3/2.3G |4.0/10.2G |13.5/14.4G |10.6/8.1G |
|CBCL|361 |2429 | 49 |0.0376 |8.9/8.8G |18.0/46.8G|45.6/49.4G |30.9/29.8G |
|CBCL|361 |2429 | 49 |0.0373 |14.6/14.5G |29.0/75.7G |84.6/91.2G |51.5/53.8G |
|ORL|10304 |400 | 25 |0.0365 |1.8/2.7G | 6.5/14.5G |9.0/9.1G|7.4/5.4G |
|ORL|10304 |400 | 25 |0.0335 |14.1/20.1G | 30.3/66.9G |98.6/9.1G|33.9/38.2G |
|ORL|10304 |400 | 25 |0.0332 |33.0/51.5G | 63.3/139.0G |256.8/193.5G|76.5/82.4G |
表での各手法の値はTime(in seconds)/FLOPSを表しており、小さいほどいいです。いずれの場合もGCD(提案手法)が高速であることがわかります。

### 実験2(疎なデータに対して)
比較する手法はGCD, FastHALS[4]、BlockPivot[6]。
もちいたデータは
1. Yahoo-News : ニュース記事、$21839\times 2340, k=20$
2. MNIST: 手書き数字の画像、$780 \times 60000, k=10$
3. RCV1: ニュース記事、$31025\times 152120, k= 15$

時間における相対誤差の推移を計測しています。
[論文7ページFigure２](http://www.cs.utexas.edu/~cjhsieh/nmf_kdd11.pdf)のグラフを見てください(重ね重ねごめんなさい…)。この場合も提案手法が高速かつ良い局所最適解に達していることがわかります。

### 実装上の注意
論文での実装はMATLABです。しかしMATLABはforループが遅いため、GCDの反復部分をmex(Cで作った関数を呼び出す)で高速化しているそうです。
公平な比較じゃない気がしますが、筆者曰くFLOPsでも比較してるから公平性はあるとのこと。あまり腑に落ちないですが。

## +α(KL-divergenceへの提案手法)
時間がないので今回は概略にととどめますが、この論文ではKL-divergenceにおける新しいアルゴリズムも提案しています。
* KL-divergence

$$
L(W, H) = \sum_{i,j}V_{i,j}log(\frac{V_{i,j}}{(WH)_{i,j}})-V_{i,j}+(WH)_{i,j}
$$

フロべニウスノルムの二乗の場合と異なり、変数選択は用いない循環座標降下法(CCD)になっています。
更新幅$s$をニュートン法を用いて高速に求めるのが特徴的な手法です。

$$
s \leftarrow \max(-W_{i,r}, s-h_{i,r}^{'}(s)/h_{i,r}^{''}(s))\\
h_{i,r}^{'}(s) = \sum_{j=1}^{n}H_{r,j}\left( 1 - \frac{V_{i,j}}{(WH)_{i,j}+sH_{r,j}} \right) \\
h_{i,r}^{''}(s) = \sum_{j=1}^{n}1 - \frac{V_{i,j}H_{r,j}^{2}}{((WH)_{i,j}+sH_{r,j})^{2}}
$$

この手法は同じ密なデータでMultiplicative Updateと比較されており、高速化に成功しています。
結果は[論文8ページTable3](http://www.cs.utexas.edu/~cjhsieh/nmf_kdd11.pdf)。
ただKL-divergenceに関しては2012年にsBCD[7]というより高速な手法が提案されています…。

## 結論
本論文では目的関数がフロべニウスノルムの二乗に用いた時に、変数選択を用いた座標降下法を提案しました。またKL-divergenceにおいてもニュートン法を利用した循環座標降下法を提案しました。この二つの手法に関して実データを用いて実験を行い、高速化に成功したことを確認できました。

## その後&感想
以上が論文の内容になります。KL-divergenceの話もそうですが、収束の保証の議論、Inner Updateの停止条件、各ステップでの計算量の議論などここに書ききれなかった話が論文にはあるので興味がある人は是非。

フロべニウスノルムの二乗を用いた際のNMFの高速化に関しては、この論文の手法(GCD)現在が最速だと思われます。そう、この場合の高速化に関してはこの論文(2011)を最後に終わってしまったのです(もしあったら是非教えて下さい！)。かっこいいですね。僕も一分野を終わらせる論文を書いてみたかったです。

ただGCDはあまり人気がなく、基本的にNMFに使用されているアルゴリズムはMultiplicative Updateが多い印象です。実装が明らかにMultiplicative Updateの方が簡単で、欲しい結果が得られているからなのかなーとか思いますが実際のところはわかりません。おしえてエロい人。

NMFの研究は現在も続けられており、毎年多くの論文が出ています。特にNMFの亜種や使用条件を絞った場合に性能を向上させる手法などが多い印象です。それには現実に欲しい結果と誤差の小ささが必ずしも関係してないことや、初期値依存性(反復更新の際にどの初期値を入れるのがベストかは基本不明な問題)など、残された問題や改良できる余地が多くあるからだと思います。とにかく今後も目が離せないですね。

来年からは少なくとも仕事で論文を読むようなことはないでしょうが、せっかく大学で論文を読む能力を身につけたのだから、年に一回こういう場でアウトプットできるぐらいには論文を読んでいきたいと思います。

## 参考文献
これを読んでNMFに興味を持った方はぜひ参考にしてみてください。
### 論文
[1]. [Fast Coordinate Descent Methods with Variable Selection for Non-negative Matrix Factorization](http://www.cs.utexas.edu/~cjhsieh/nmf_kdd11.pdf) ：今回紹介した論文です。コードは[筆者のホームページ](http://www.cs.utexas.edu/~cjhsieh/nmf/)に公開されています。
[2]. [Learning the parts of objects by non-negative matrix factorization](https://www.nature.com/articles/44565):NMFが最初に世に出た時の論文です。
[3]. [Algorithms for Non-negative Matrix Factorization](http://papers.nips.cc/paper/1861-algorithms-for-non-negative-matrix-factorization.pdf)：一番有名なNMFのアルゴリズムであるMultiplicative Updateの論文です。
[4].[Fast local algorithms for large scale nonnegative matrix and tensor factorizations](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.214.6398&rep=rep1&type=pdf):FastHALSの論文です。
[5]. [Projected gradient methods for nonnegative matrix factorization](https://www.mitpressjournals.org/action/captchaChallenge?redirectUrl=https%3A%2F%2Fwww.mitpressjournals.org%2Fdoi%2Fpdfplus%2F10.1162%2Fneco.2007.19.10.2756%3Fcasa_token%3DjlubOQwKdb4AAAAA%253AHBNw9iRU6Xsf2grl1CRBhiopkx7uHNNp1b5mnQHlrQ38fgbqZU-BZPoyoLerOTeUtEfc7_GqDC9_Pw)：実験の比較に出てきた手法です。
[6]. [Toward Faster Nonnegative Matrix Factorization: A New Algorithm and Comparisons](https://smartech.gatech.edu/bitstream/handle/1853/25538/GT-CSE-08-03.pdf?sequence=1&isAllowed=y)：実験の(ry
[7]. [Fast Bregman Divergence NMF using Taylor Expansion and Coordinate Descent](https://www.ime.usp.br/~jstern/miscellanea/seminario/Lebanon13.pdf):2012に出てきた手法です。僕はこの手法はFastHALSの一般化と解釈しています。

### 日本語記事
[非負値行列因子分解](http://www.kecl.ntt.co.jp/people/kameoka.hirokazu/publications/Kameoka2012SICE09published.pdf)：NTT研究所の方によるNMFの解説記事です。非常にわかりやすいです。
[Nonnegative Matrix Factorization:NMFの最適化法についてのまとめ](http://kkimura.hatenablog.com/entry/2017/10/07/211926)：色々なNMFのアルゴリズムの概説が載っています。

誤字脱字の指摘から記事の内容に関するツッコミなどなんでも歓迎です！
twitterアカウント[localmin](https://twitter.com/localmin?lang=ja)やメールアドレス localmin9201009"at"gmail.comまでお願いします！