諸々

信号処理・機械学習に興味を持っているエンジニア見習いのメモです。技術的なことから趣味のことを書いていきます。

非負値行列因子分解(NMF : Non-negative Matrix Factorization)について

はじめに

以前Qiitaに投稿したものをメモ用にとこちらに移植しました.勉強ついでに少しイジって更新しました. Support Vector Machine (SVM)や最近ではDeep Neural Network (DNN) などが流行っていますが,「こんな手法もありますよ」というのを知ってほしくて書きました.

非負値行列因子分解(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参照)

Itakura–Saito distance - Wikipedia

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

他にもこれら3つの距離関数を表現している,βダイバージェンスなどもあったりします. これらの距離アルゴリズムを乗法的更新アルゴリズムとして使用します. 細かい式の導出や変形は省略します.

変形後の更新式は,

ユークリッド距離(EUC

 \displaystyle
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)

 \displaystyle
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)

 \displaystyle
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 MIT Press Journals

実装

実際にPythonで実装した例を提示します.

開発環境

Python 2.7.11(Qiitaに挙げてから更新してない・・・)

numpy 1.11.1

cmake 3.6.2

# 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仕様のものだけを記載しています.

テスト

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

#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++版はこちら

おわりに

Qiitaに挙げてから随分時間が立って少し忘れていたので,よい復習になりました.

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