非負矩陣分解（NMF）原理及算法實現

一、矩陣分解回想

矩陣分解是指將一個矩陣分解成兩個或者多個矩陣的乘積。對於上述的用戶-商品（評分矩陣），記為能夠將其分解為兩個或者多個矩陣的乘積，如果分解成兩個矩陣和 。我們要使得矩陣和 的乘積能夠還原原始的矩陣

當中，矩陣表示的是m個用戶於k個主題之間的關系，而矩陣表示的是k個主題與n個商品之間的關系

通常在用戶對商品進行打分的過程中，打分是非負的，這就要求：

這便是非負矩陣分解（NMF）的來源。

二、非負矩陣分解

2.1、非負矩陣分解的形式化定義

上面介紹了非負矩陣分解的基本含義。簡單來講，非負矩陣分解是在矩陣分解的基礎上對分解完畢的矩陣加上非負的限制條件。即對於用戶-商品矩陣找到兩個矩陣和 ，使得：

同一時候要求：

2.2、損失函數

為了能夠定量的比較矩陣和的近似程度，提出了兩種損失函數的定義方式：
歐幾里得距離：

KL散度：

在KL散度的定義中，。當且僅當時取得等號。

當定義好損失函數后，須要求解的問題就變成了例如以下的形式，相應於不同的損失函數：

求解例如以下的最小化問題：

2.3、優化問題的求解

乘法更新規則，詳細操作例如以下：

對於歐幾里得距離的損失函數：

對於KL散度的損失函數：

上述的乘法規則主要是為了在計算的過程中保證非負，而基於梯度下降的方法中，加減運算無法保證非負。事實上上述的懲罰更新規則與梯度下降的算法是等價的。以下以平方距離為損失函數說明上述過程的等價性：

平方損失函數能夠寫成：

使用損失函數對求偏導數：

依照梯度下降法的思路：

即為：

令，即能夠得到上述的乘法更新規則的形式。

2.4、非負矩陣分解的實現

from numpy import *
from pylab import *
from numpy import *

def load_data(file_path):
    f = open(file_path)
    V = []
    for line in f.readlines():
        lines = line.strip().split("\t")
        data = []
        for x in lines:
            data.append(float(x))
        V.append(data)
    return mat(V)

def train(V, r, k, e):
    m, n = shape(V)
    #先隨機給定一個W、H，保證矩陣的大小
    W = mat(random.random((m, r)))
    H = mat(random.random((r, n)))
#K為迭代次數
    for x in range(k):
        #error
        V_pre = W * H
        E = V - V_pre
        #print E
        err = 0.0
        for i in range(m):
            for j in range(n):
                err += E[i,j] * E[i,j]
        print(err)
        data.append(err)

        if err < e:
            break
#權值更新
        a = W.T * V
        b = W.T * W * H
        #c = V * H.T
        #d = W * H * H.T
        for i_1 in range(r):
            for j_1 in range(n):
                if b[i_1,j_1] != 0:
                    H[i_1,j_1] = H[i_1,j_1] * a[i_1,j_1] / b[i_1,j_1]

        c = V * H.T
        d = W * H * H.T
        for i_2 in range(m):
            for j_2 in range(r):
                if d[i_2, j_2] != 0:
                    W[i_2,j_2] = W[i_2,j_2] * c[i_2,j_2] / d[i_2, j_2]

    return W,H,data




if __name__ == "__main__":
    #file_path = "./data_nmf"
    # file_path = "./data1"
    data = []
    # V = load_data(file_path)
    V=[[5,3,2,1],[4,2,2,1,],[1,1,2,5],[1,2,2,4],[2,1,5,4]]
    W, H ,error= train(V, 2, 100, 1e-5 )
    print (V)
    print (W)
    print (H)
    print (W * H)
    n = len(error)
    x = range(n)
    plot(x, error, color='r', linewidth=3)
    plt.title('Convergence curve')
    plt.xlabel('generation')
    plt.ylabel('loss')
    show()

 這里需要注意訓練時r值的選擇：r可以表示和主題數或者你想要的到的特征數

K值的選擇：k表示訓練的次數，設置的越大模型的擬合效果越好，但是具體設置多少，要根據性價比看，看誤差曲線的變化