一年ほど前に理解した非負値行列因子分解（NMF）を再実装して確認してみようとした記事です．
Support Vector Machine(SVM)や最近ではDeep Learningなどが流行っていますが，こんな手法もありますよ，といった紹介になればと・・・
C++とPython2.7で久しぶりに実装してみました．

【2020/11/10 更新】

Python3系にアップデートしました。
Python版NMFコードをリファクタリングしました。

【2020/02/05 更新】

NMFを使ったリアルタイム音源分離のモックを作成しました。
こちら

【2018/03/23 更新】

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

開発環境

Mac OS X 10.11.6 Elcapitan

Python 3.7

numpy 1.11.1

cmake 3.6.2

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

NMFは，非負値(>=0)の元行列 V を他の２つの非負値な行列 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

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

変形後の更新式は，

ユークリッド距離（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


import numpy as np
from typing import List, Tuple


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

    def setParam(self, k: int, row: int, column: int):
        """NMFのパラメータ設定
        Args:
            k (int): 因子数
            row (int): 列数
            column (int): 行数
        """
        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: int, data: List):
        """辞書行列へのデータ設定
        Args:
            index (int): 因子インデックス ( 0 <= index < k)
            data (List): データ
        """
        if index >= self.__k and len(data) != self.__row:
            print("Please NMF.setParam(k,row,column)")
            print(f"k = {self.__k}")
            print(f"row = {self.__row}")
            return

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

    def setAnalyzData(self, data: List, k: int):
        """分解対象行列データを登録
        Args:
            data (List): 分解対象行列
            k (int): 因子数
        """
        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: int = 200) -> Tuple[np.array, np.array]:
        """テンプレートありのEUC-divergence仕様の分離処理
        Args:
            iter (int, optional): 反復更新回数. Defaults to 200.
        Returns:
            Tuple[np.array, np.array]: [辞書行列, 励起行列]
        """
        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: int = 200) -> Tuple[np.array, np.array]:
        """テンプレートなしのEUC-divergence仕様の分離処理
        Args:
            iter (int, optional): 反復更新回数. Defaults to 200.
        Returns:
            Tuple[np.array, np.array]: [辞書行列, 励起行列]
        """
        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

from NMF import NMF


"""
テスト用のスクリプト
使い方記載
"""
if __name__ == "__main__":

    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_kl_with_template(iter=200)
    # dic,act = nmf.separate_is_with_template(iter=200)

    # dic,act = nmf.separate_euc_without_template(iter=200)
    # dic,act = nmf.separate_kl_without_template(iter=200)
    # dic,act = nmf.separate_is_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などが注目されていますが，こんな簡単なアルゴリズムもあるんだ！っという紹介になればいいと思ってます．

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