Edited at

非負値行列因子分解(NMF : Non-negative Matrix Factorization)を改めてやり直してみた

More than 1 year has passed since last update.

一年ほど前に理解した非負値行列因子分解(NMF)を再実装して確認してみようとした記事です.

Support Vector Machine(SVM)や最近ではDeep Learningなどが流行っていますが,こんな手法もありますよ,といった紹介になればと・・・

C++とPython2.7で久しぶりに実装してみました.

【2018/03/23 更新】

一部本文を修正しました. 加えてJuliaでの実装も行いましたのでそちらのリンクを貼りました.


開発環境

Mac OS X 10.11.6 Elcapitan

Python 2.7.11(そろそろ3系に上げないと・・・)

numpy 1.11.1

cmake 3.6.2


非負値行列因子分解(Non-negative Matrix Factorization)

NMFは,非負値(>=0)の元行列 V を他の2つの非負値な行列 W, Hの積で近似するアルゴリズムです.

このアルゴリズム自体は,特徴抽出や自然言語処理,クラスタリング手法として用いられます.

またその非負値性から,音響信号処理におけるパワースペクトルを用いたものもあり,スペクトログラムの中から特定の楽器の音色の発音時刻を見つけるといった使い方もされています.


数式表現

i * jのもと行列に対しての,NMFの数式は,


k \in R \\

V \approx WH\\

V \in R^{i \times j}\\

W \in R^{i \times k} \qquad H \in R^{k \times j}

となります.

入力されてくる行列Vの中身は分かっていますが,WとHの中身は不明なので,NMFではこのVとWHの距離を最小化します.


更新アルゴリズム

NMFでよく知られている(そして簡単)なアルゴリズムは,乗法的更新アルゴリズムです.

上でも述べましたが,XとWHの距離(誤差)を定義し,それを最小化しなければなりません.

NMFで使われる距離は私が知っているもので3種類存在します.


  • ユークリッド距離(EUC)
    よく使われる距離関数.

D_{EUC}(V, WH) = (V-WH)^2


  • 一般化カルバック・ライブラー距離(KL)
    2つの確率分布の違いをはかる尺度だったりします."距離"というわけではないらしいですが,KLダイバージェンスを最小化することと,尤度の最大化を同じ意味らしいです.

D_{KL}(V, WH) = V log \frac{V}{WH} - V + WH


  • Itakura-Saito距離(IS)
    音響処理におけるスペクトルの近似の差を示す尺度らしいです.
    知覚的尺度ではないが,知覚的類似性を反映することを意図してるらしい・・・.(wiki参照)



https://en.wikipedia.org/wiki/Itakura–Saito_distance



D_{IS}(V,WH) = \frac{V}{WH} - log \frac{V}{WH} -1

他にもこれら3つの距離関数を表現している,βダイバージェンスなどもあったりします.

これらのアルゴリズムを乗法的更新アルゴリズムとして使用します.

細かい式の導出や変形は省略します.

変形後の更新式は,

* ユークリッド距離(EUC)


H_{n+1} \gets H_{n} \frac{W^{T}_{n}V_{n}}{W^{T}_{n}W_{n}H_{n}}, \qquad W_{n+1} \gets W_{n} \frac{V_{n}H^{T}_{n+1}}{W_{n}H_{n+1}H^{T}_{n+1}}


  • 一般化カルバック・ライブラー距離(KL)


H_{n+1} \gets H_{n} \frac{W^{T}_{n}\frac{V_{n}}{W_{n}H_{n}}}{W^{T}_{n}}, \qquad W_{n+1} \gets W_{n} \frac{\frac{V_{n}}{W_{n}H_{n+1}}H^{T}_{n+1}}{H^{T}_{n+1}}


  • Itakura-Saito距離(IS)


H_{n+1} \gets H_{n} \sqrt{\frac{\frac{W_{n}}{W_{n}H_{n}}^{T}\frac{V_{n}}{W_{n}H_{n}}}{\frac{W_{n}}{W_{n}H_{n}}^{T}}}, \qquad W_{n+1} \gets W_{n} \sqrt{\frac{\frac{V_{n}}{W_{n}H_{n+1}}\frac{H_{n+1}}{W_{n}H_{n+1}}^{T}}{\frac{H_{n+1}}{W_{n}H_{n+1}}^{T}}}

こんな感じだったかと・・・(若干うろ覚え.間違ってたら申し訳ない)

気をつけないといけないのは,数式の中に入ってるHが{n+1}になっていること!

これを忘れると,全く収束しない事態に陥ります.

これを反復的に実行することで,VとWHの距離を小さくしていきます..

参考資料・URL



http://papers.nips.cc/paper/1861-algorithms-for-non-negative-matrix-factorization.pdf

http://www.mitpressjournals.org/doi/pdf/10.1162/neco.2008.04-08-771




実装


NMF.py


# coding:utf-8

import numpy as np

class NMF():
"""
簡単なNMFを行うクラス
"""

def setParam(self, k, row, column):
"""
:param k
:param row: 列
:param column: 行
:return:
"""

self.__k = k
self.__row = row
self.__column = column

self.__dictionary = np.random.random_sample([self.__row, self.__k])
self.__activation = np.random.random_sample([self.__k, self.__column])

def setDictionary(self, index, data):
"""
テンプレートありの場合,テンプレートをセットしてください(デフォルトは乱数)
:param index: どのテンプレートかのindex (0 <= index < k)
:param data: テンプレートデータ
:return:
"""

if index >= self.__k and len(data) != self.__row:
print "Please NMF.setParam(k,row,column)"
print "k = " + str(self.__k)
print "row = " + str(self.__row)
return

self.__dictionary[:, index] = np.array(data[:self.__row], np.float32)

def setAnalyzData(self, data, k):
"""
一番最初に,このメソッドを呼んでください
:param data: 分解するデータ
:param k
:return:
"""

if len(np.shape(data)) == 1:
self.__data = np.ones([np.shape(data)[0], 1], np.float32)
self.setParam(k, np.shape(data)[0], 1)
else:
self.__data = data
self.setParam(k, np.shape(data)[0], np.shape(data)[1])

def separate_euc_with_template(self,iter=200):
"""
テンプレートありのEUC仕様の分離処理を行う
:param iter:
:return:
"""

counter = 0

while counter < iter:
approx = np.dot(self.__dictionary , self.__activation)

wh = np.dot(np.transpose(self.__dictionary) , self.__data)
wt = np.dot(np.transpose(self.__dictionary) , approx)

bias = wh/wt
bias[np.isnan(bias)] = 0

self.__activation = self.__activation * bias
counter += 1
return self.__dictionary,self.__activation

def separate_euc_without_template(self,iter=200):
"""
テンプレートなしのEUC仕様の分離処理を行う
:param iter:
:return:
"""

counter = 0

while counter < iter:
approx = np.dot(self.__dictionary , self.__activation)

wh = np.dot(np.transpose(self.__dictionary) , self.__data)
wt = np.dot(np.transpose(self.__dictionary) , approx)

bias = wh/wt
bias[np.isnan(bias)] = 0

self.__activation = self.__activation * bias

approx = np.dot(self.__dictionary,self.__activation)
wh = np.dot(self.__data,np.transpose(self.__activation))
wt = np.dot(approx,np.transpose(self.__activation))

bias = wh/wt
bias[np.isnan(bias)] = 0
self.__dictionary = self.__dictionary * bias
counter += 1

return self.__dictionary,self.__activation


NMF自体を実行するクラスはこんな感じで今回実装しました.

文量の問題で今回は,EUC仕様のものだけを記載しています.

なかなか面倒なコードになっていますが,気にしないで頂ければ・・・.


テスト

テスト用のコード(使い方含む)も作ったので実行してみました.


main.py


#coding:utf-8

import numpy as np

"""
面倒くさいけれど、この2つをimportしないといけない
"""

from NMF import NMF

"""
テスト用のスクリプト
使い方記載
"""

if __name__ == "__main__":

"""
基本的にはコンストラクタで面倒なk,row,columnを渡さなくてもOK
適当に作っただけなので・・・
"""

nmf = NMF()

"""
コンストラクタを作ったあとには、まずこれを呼んでください!
ここで、k,row,columnの設定をするので、これを呼ばないとsetDictionaryが動かなくなります!
"""

nmf.setAnalyzData([[1,2,3,4],[2,3,4,5]],k=3)

"""
ここでテンプレートをセットする
"""

nmf.setDictionary(0,[0.0,2.0])
nmf.setDictionary(1,[1.0,6.0])
nmf.setDictionary(2,[11.0,10.0])

"""
NMF開始
引数には、アルゴリズムと反復更新回数を渡しておく
"""

dic,act = nmf.separate_euc_with_template(iter=200)
# dic,act = nmf.separate_euc_without_template(iter=200)
"""
結果表示
"""

print "=========================== Dictionary ==========================="
print dic
print "=========================== Activation ==========================="
print act

"""
ちゃんと分解できているかの確認
"""

print "=========================== Approx ==========================="
print np.dot(dic,act)


このテストでは,


V = \begin{pmatrix}
1 & 2 & 3 & 4\\
2 & 3 & 4 & 5
\end{pmatrix}
\qquad k=3

といった条件の行列に対してNMFを行っている.

この行列に対して,setDictionary()でセットした辞書行列


W = \begin{pmatrix}
0 & 1 & 11 \\
2 & 6 & 10
\end{pmatrix}

を使って,励起行列を更新します.


テスト結果

main.pyスクリプトを実行すると・・・


$ python main.py
=========================== Dictionary ===========================
[[ 0. 1. 11.]
[ 2. 6. 10.]]
=========================== Activation ===========================
[[ 0.11776982 0.05902263 0.00976704 0.33375605]
[ 0.168019 0.20895539 0.24616295 0.1367387 ]
[ 0.07563464 0.16282224 0.25034882 0.35120557]]
=========================== Approx ===========================
[[ 1. 2. 3. 4.]
[ 2. 3. 4. 5.]]

ちゃんと,励起行列は更新されてVとWHの距離が0になりました.


ソースコード

今回は省略しました,KL・IS版も含めたものをGitHub上で公開しています.

また,今回は紹介を省略したC++版のソースコード(構成はPython版とほぼ同じ)も同様に公開しています.

Python版はこちら

C++版はこちら

Julia版はこちら


最後に・・・

約一年前に勉強したことを久しぶりに文章・数式にするのはよい復習になりました.

はじめにでも書きましたが,最近ではDNNなどが注目されていますが,こんな簡単なアルゴリズムもあるんだ!っという紹介になればいいと思ってます.

今度は,これを使った実践的な例を紹介しようかなぁと考えています.