NMF非负矩阵分解

    NMF非负矩阵分解是将一个非负矩阵分解成两个非负矩阵,处理有些实际问题时数据往往是非负数,其他的矩阵分解如SVD分解有时会将矩阵分解成含有负数矩阵,负数出现与实际情况相悖。NMF在图像处理、推荐系统、数据降维中有着广泛的应用,NMF在实现方法上也与其他矩阵分解有质的区别。

一、理解NMF

   NMF将一个非负矩阵V分解成两个非负矩阵W、V,三者的关系是:

V≈W*H

注意中间是约等于的关系,NMF需要找出两个非负矩阵来近似原来的矩阵。不妨设V是m行n列的矩阵:V∈Rm*n,W是m行r列矩阵:WRm*r,H是r行n列矩阵:HRr*n。从几何角度看,V代表了一个n维线性空间,W代表了一个r维线性空间,由于r<n,所以W是V的子空间,H矩阵代表系数矩阵,也可以理解为坐标矩阵,所以也可以这样定义NMF:找到一个线性空间V的子空间W,使得空间V中的向量投影到子空间W中后误差最小,如下图:

几何意义.png

v1,v2,v3是V空间的向量,e1,e2,e3v1,v2,v3投影到W子空间后的误差,满足所有误差最小子空间W即为目标矩阵,在NMF算法中W称为基矩阵,H称为系数矩阵。

    来看一个具体的例子,下面的程序根据10位受访者对10部影片打分,利用NMF可分别对用户和电影进行分类:

import numpy as np
from sklearn.decomposition import NMF
from _collections import defaultdict
import  matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
if __name__=='__main__':
    ratings = [
        [5.0, 5.0, 4.5, 4.5, 5.0, 3.0, 2.0, 2.0, 0.0, 0.0],
        [4.2, 4.7, 5.0, 3.7, 3.5, 0.0, 2.7, 2.0, 1.9, 0.0],
        [2.5, 0.0, 3.3, 3.4, 2.2, 4.6, 4.0, 4.7, 4.2, 3.6],
        [3.8, 4.1, 4.6, 4.5, 4.7, 2.2, 3.5, 3.0, 2.2, 0.0],
        [2.1, 2.6, 0.0, 2.1, 0.0, 3.8, 4.8, 4.1, 4.3, 4.7],
        [4.7, 4.5, 0.0, 4.4, 4.1, 3.5, 3.1, 3.4, 3.1, 2.5],
        [2.8, 2.4, 2.1, 3.3, 3.4, 3.8, 4.4, 4.9, 4.0, 4.3],
        [4.5, 4.7, 4.7, 4.5, 4.9, 0.0, 2.9, 2.9, 2.5, 2.1],
        [0.0, 3.3, 2.9, 3.6, 3.1, 4.0, 4.2, 0.0, 4.5, 4.6],
        [4.1, 3.6, 3.7, 4.6, 4.0, 2.6, 1.9, 3.0, 3.6, 0.0],
    ]
    movie={
        1:'星球大战',
        2:'黑客帝国',
        3:'盗梦空间',
        4:'哈利波特',
        5:'霍比特人',
        6:'六壮士',
        7:'拯救大兵瑞恩',
        8: '兵临城下',
        9: '魔窟寻谍',
        10: '大逃亡'
    }
    A=np.asmatrix(ratings,dtype=np.float)
    max_components=2
    nmf=NMF(n_components=max_components, init='random', random_state=1)
    A_dash=nmf.fit_transform(A)
    F=nmf.components_
    xf,yf=F[0,:],F[1,:]
    for i in range(A_dash.shape[0]):
        print('用户id:%d 成分1得分:%0.2f 成分1得分:%0.2f'%(i+1,A_dash[i,0],A_dash[i,1]))
        print()
    print('='*50)
    for i in range(F.shape[1]):
        print('%s 成分1得分:%0.2f 成分1得分:%0.2f'%(movie[i+1],F[0,i] ,F[1,i] ))
        print()
    plt.figure(1)
    plt.title('基矩阵')
    x,y=A_dash[:,0],A_dash[:,1]
    plt.scatter(x,y)
    plt.xlabel('成分1得分')
    plt.ylabel('成分2得分')
    for i in range(A_dash.shape[0]):
        plt.annotate(i+1, (A_dash[i, 0], A_dash[i, 1]) )
    plt.figure(2)
    plt.title('系数矩阵')
    plt.scatter(xf, yf)
    plt.xlabel('成分1得分')
    plt.ylabel('成分2得分')
    for i in range(F.shape[1]):
        plt.annotate(movie[i+1],(F[0,i] ,F[1,i] ))
    plt.show()

ratings变量代表受访者的打分信息,是一个10*10的非负矩阵,上述代码利用nmf将ratings分解成两个矩阵,一个是矩阵A_dash∈R10*2,A_dash对应NMF分解的基矩阵,打印A_dash有以下信息:

w得分.png

根据A_dash可以对用户进行分类,比如id=1的用户,根据打分的情况可以判断出他偏好成分2的电影,但不喜欢成分1的电影,再比如id=3的用户偏好成分1的电影,而成分1和成分2分别对应哪些题材电影呢?这可以通过系数矩阵收集到相应信息,系数矩阵对应代码中矩阵F,打印F矩阵信息如下:

H得分.png

可以看出星球大战成分2得分较多,id1和id2的用户都喜欢星球大战、黑客帝国、霍比特人这类电影,这当然不是偶然的,可以发现成分2得分较多的电影大都是科幻、魔幻类电影,而成分1得分较多的是战争题材电影。这样一来就通过NMF对受访者以及电影实现了分类,如下图:

Figure_1.pngFigure_2.png

还可以利用NMF开发推荐系统,比如一个用户没有对其中一个偏好成分1(战争题材电影)的电影打分,但通过用户对其他电影打分可以判断出该用户偏向于成分1的电影,那么就可以向其推荐战争题材的电影。

二、NMF公式推导

2.1、矩阵范数

    分析向量、矩阵时,常常利用范数将向量、矩阵变为一个实数,变为实数后就可以利用多种手段如微分、积分、幂级数等手段做进一步处理。设有一个矩阵A∈Rm*n,矩阵A的任意一个元素用aij表示,常用矩阵范数形式有1-范数:

范数一.gif

Frobenius范数或称F-范数:

f范数.gif

范数几何意义是向量或矩阵的长度,当线性空间引入范数后也就有了距离的概念,这很好理解,比如两个向量v和u,v-u代表两者之间一个向量,而范数||v-u||就表示v和u之间向量长度,也就得到了v和u之间距离,所以线性赋范空间是一个距离空间。距离空间与熟识欧几里得空间还差一个角度的概念,角度需要在距离空间的基础上再引入内积。

    基矩阵W与系数矩阵H相乘后的矩阵用B表示,V-B表示分解后矩阵B与原矩阵V的误差,结合上面的分析,范数||V-B||较小时说明B与原矩阵有较好的相似性,NMF算法的目标可以描述为求下面F-范数平方的最小值:

nmf范数.gif              

W是m行r列矩阵:WRm*r,H是r行n列矩阵:HRr*n,则bij可表示为:

bij.gif                                           

②代入①式后可得下面的损失函数:

nmf范数2.gif

公式(3)中共有m*n个未知参数,可以通过梯度下降法求未知参数,上式中由于含有平方项,为后期处理方便通常把损失函数表示成下面的形式:

范数2.png

2.2 、最大似然法推导损失函数

    上面利用范数导出导出了损失函数(3.1)式,接下来利用最大似然法从另一个角度理解NMF算法,依然利用之前的假设:基矩阵W与系数矩阵H相乘后的矩阵用B表示,则V=B+E,E代表了V与B之间的误差,具体到矩阵每个元素可以这样表示:

eij.gif

假设eij为数学期望为0,方差为σ2的正态分布,即正态分布.gif,最大似然法即为求下面函数最大值:

最大似然法.png

由于方差σ大于0,所以求上式最大值实质为求最狭窄.gif最小值,这与用矩阵范数导出的损失函数是一致的,这个结论并不是偶然的,矩阵法本质是最小二乘,而最小二乘法、线性回归都是将误差设定为正态分布。NMF在应用中根据实际情况,误差项也可以用其他概率分布实现,之前文章曾介绍过,最大似然法与求交叉熵最小值是本质一个问题,如果将泊松分布与熵模型结合可以得到NMF另一种损失函数形式-KL散度损失函数,形式如下:

泊松分布.png

本篇重点研究矩阵分解后误差为正态分布的情况,即NMF损失函数为公式(3.1),接下来对式(3.1)求导得到各个参数梯度。

2.3、求损失函数最优解

    设wit为基矩阵W(W∈Rm*r,t≤r)任意一个元素,由链式法则有:

链式法则.gif

公式②可得:

htj.gif

而当确定w的下标i和t后,观察可以发现:wit只与矩阵B=W*H的第i行向量有关,所以有:

求导1.png

这样得到未知参数wit的梯度为:

wt梯度.gif    

注意,④式中i和t是常量,j,k是迭代变量,通常式写成下面矩阵相乘的形式:

梯度1.gif        (4.1)

同理,可以得到系数矩阵H任意一个元素htj的梯度:

梯度2.gif       (4.2)

求目标函数最小值,则需要按梯度反方向更新参数:

推导公式.png

上式中μit 、 ηtj称为步长系数,在神经网络中也成为学习率,可以发现NMF算法中每个未知参数的步长系数都是不同的,最重要的一点,NMF必须确保分解后矩阵是非负数,μit 、 ηtj取数是受到一定限制的,可以通过下面的换元法保证分解后矩阵的非负性:

步长西湖水.png

可得到未知参数最终的迭代计算公式:

迭代公式43.png

通过换元法变换步长系数不仅能确保矩阵分解后非负,还能简化计算变为矩阵乘法。

三、NMF代码实现

    利用以上推导过程实现NMF,示例代码通过公式(4.3)求出(3.1)最优解后得到未知参数,并将分解后矩阵与sklearn自带的NMF算法比较,通过矩阵范数计算可以发现两者准确率是相同的。

-请登陆后阅读余下文章-
登录|注册