推薦算法——非負矩陣分解(NMF)

一、矩陣分解回想

在博文推薦算法——基於矩陣分解的推薦算法中，提到了將用戶-商品矩陣進行分解。從而實現對未打分項進行打分。

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

Vm×n≈Wm×k×Hk×n=V^m×n

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

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

Wm×k⩾0

Hk×n⩾0

這便是非負矩陣分解(Non-negtive Matrix Factorization, NMF)的來源。

二、非負矩陣分解

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

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

Vm×n≈Wm×k×Hk×n=V^m×n

同一時候要求：

Wm×k⩾0

Hk×n⩾0

2.2、損失函數

為了能夠定量的比較矩陣 Vm×n 和矩陣 V^m×n 的近似程度。在參考文獻1中作者提出了兩種損失函數的定義方式：

• 平方距離

∥A−B∥2=∑i,j(Ai,j−Bi,j)2

• KL散度

D(A∥B)=∑i,j(Ai,jlogAi,jBi,j−Ai,j+Bi,j)

在KL散度的定義中， D(A∥B)⩾0 。當且僅當 A=B 時取得等號。

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

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

• minimize∥V−WH∥2s.t.W⩾0,H⩾0

• minimizeD(V∥WH)s.t.W⩾0,H⩾0

2.3、優化問題的求解

在參考文獻1中，作者提出了乘法更新規則(multiplicative update rules)，詳細的操作例如以下：

對於平方距離的損失函數：

Wi,k=Wi,k(VHT)i,k(WHHT)i,k

Hk,j=Hk,j(WTV)k,j(WTWH)k,j

對於KL散度的損失函數：

Wi,k=Wi,k∑uHk,uVi,u/(WH)i,u∑vHk,v

Hk,j=Hk,j∑uWu,kVu,j/(WH)u,j)∑vWv,k

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

平方損失函數能夠寫成：

l=∑i=1m∑j=1n[Vi,j−(∑k=1rWi,k⋅Hk,j)]2

使用損失函數對 Hk,j 求偏導數：

∂l∂Hk,j=∑i=1m∑j=1n[2(Vi,j−(∑k=1rWi,k⋅Hk,j))⋅(−Wi,k)]=−2[(WTV)k,j−(WTWH)k,j]

則依照梯度下降法的思路：

Hk,j=Hk,j−ηk,j∂l∂Hk,j

即為：

Hk,j=Hk,j+ηk,j[(WTV)k,j−(WTWH)k,j]

令 ηk,j=Hk,j(WTWH)k,j ，即能夠得到上述的乘法更新規則的形式。

2.4、非負矩陣分解的實現

對於例如以下的矩陣：

通過非負矩陣分解。得到例如以下的兩個矩陣：

對原始矩陣的還原為：

實現的代碼

#!/bin/python

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 = mat(random.random((m, r)))
    H = mat(random.random((r, n)))

    for x in xrange(k):
        #error 
        V_pre = W * H
        E = V - V_pre
        #print E
        err = 0.0
        for i in xrange(m):
            for j in xrange(n):
                err += E[i,j] * E[i,j]
        print 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 xrange(r):
            for j_1 in xrange(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 xrange(m):
            for j_2 in xrange(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 


if __name__ == "__main__":
    #file_path = "./data_nmf"
    file_path = "./data1"

    V = load_data(file_path)
    W, H = train(V, 2, 100, 1e-5 )

    print V
    print W
    print H
    print W * H

收斂曲線例如以下圖所看到的：

''' Date:20160411 @author: zhaozhiyong '''

from pylab import *
from numpy import *

data = []

f = open("result_nmf")
for line in f.readlines():
    lines = line.strip()
    data.append(lines)

n = len(data)
x = range(n)
plot(x, data, color='r',linewidth=3)
plt.title('Convergence curve')
plt.xlabel('generation')
plt.ylabel('loss')
show()

參考文獻

• Algorithm for Non-negative Matrix Factorization

• 白話NMF（Non-negative Matrix Factorization）——Matlab 實現