集体智慧编程学习之非负矩阵因式分解

欢迎关注我的个人博客blog.timene.com

大学时我的线性代数老师寿继麟，当时六十多岁带着一副金丝眼镜精神矍铄，传说是我最尊敬的余德鴻副校长的老师。上课的课本是寿老师写的打印出来给大家，很便宜。我虽然不好好学习，但是在好老师的带领下，也不至于学的太差。余校长隔段时间总会给我们上一堂课，虽然没有说过我们什么，但总能让我们惭愧进而刻苦上一段时间。记得他说过这样一个问题，会场人声嘈杂，给你两个麦克风，怎么在一无所知的情况下分离出你想听到的那个人的声音。十年后的今年我才真正开始考虑这个问题，余校长提出的问题，原来寿老师当时就帮我们解决了。

非负矩阵因式分解，非负是说矩阵的元素都是正数或者零，因式分解就是对矩阵A，寻找矩阵B和C，使得A=B*C。这样分解有什么用呢？还是为了解决分类的问题。区别人声，也是一个分类问题。分类的问题在之前《集体智慧编程学习之分类系统》中提到了贝叶斯和费舍尔，在《集体智慧编程学习之聚类系统》中提到了K-均值聚类，贝叶斯和费舍尔需要事前定义好分类，然后训练模型，K-均值聚类的计算量还是蛮大的。而非负矩阵因式分解可以一次把所有的相关性计算出来，还是很有必要好好学习的。

矩阵因式分解的目的呢，就是把矩阵B拆分成为两个矩阵相乘，比如矩阵F和矩阵W，这样B=W*F。矩阵的乘法当时寿老师教的很清楚，第一个矩阵（比如W）的列必须于第二个矩阵（比如F）的行数相等才能相乘，乘积矩阵（比如B）的每个元素的取值，是通过将矩阵W中相同行的值与矩阵F中相同列的值相乘，然后将乘积相加所得。（这里我不举例子了，随便Google就有了）。这里将B拆分为F和W的乘积，F代表特征矩阵，W代表权值矩阵，下面详细说说这两个东西。

F特征矩阵：在该矩阵中，每一个特征对应一行，每个单词对应一列，矩阵中的数字代表了某个单词相对于某个特征的重要程度；

W权值矩阵：在改矩阵中，每一行对应一本书，每一列对应一个特征，矩阵中的数字代表了将某个特征应用于某篇文章的程度；

这样W*F就可以重构出B，一行对应一本书，一列对应一个单词。

通过计算最佳的特征矩阵和权值矩阵，算法尝试尽可能大地来重新构造文章矩阵。这里我们先定义一种方法来衡量最终结果与理想结果的接近程度：

def difcost(a,b):  dif=0  for i in range(shape(a)[0]):    for j in range(shape(a)[1]):      # Euclidean Distance      dif+=pow(a[i,j]-b[i,j],2)  return dif

函数difcost针对两个同样大小的矩阵遍历其中的每个值，并将两者的间差值的平方累加起来。

下面我们需要一种方法能够逐渐地更新矩阵，以使成本函数的计算值逐步降低。《集体智慧编程学习之优化系统》里的退火法和遗传算法都能满足要求。这里我来学习一种更为的方法，乘法更新矩阵法则。法则的具体原理先不学习，先直接拿来用。有兴趣的看这篇论文《Algorithms for Non-negative Matrix Factorization》

这个法则产生四个新的更新矩阵：

hn：转置后的权重矩阵与数据矩阵相乘而来；

hd：转置后的权重矩阵与原权重矩阵相乘，再与特征矩阵相乘而来；

wn：数据矩阵与转置后的特征矩阵相乘而来；

wd：权重矩阵与特征矩阵相乘，再与转置后的特征矩阵相乘而来；

更新特征矩阵和权重矩阵，我们将特征矩阵中的每个值与hn中的对应值相乘，并除以hd中的对应值；再将权重矩阵中的每个值与wn中的对应值相乘，并除以wd中的对应值下面看用python的算法实现：

def factorize(v,pc=10,iter=50):  ic=shape(v)[0]  fc=shape(v)[1]  # Initialize the weight and feature matrices with random values  w=matrix([[random.random() for j in range(pc)] for i in range(ic)])  h=matrix([[random.random() for i in range(fc)] for i in range(pc)])  # Perform operation a maximum of iter times  for i in range(iter):    wh=w*h        # Calculate the current difference    cost=difcost(v,wh)        if i%10==0: print cost        # Terminate if the matrix has been fully factorized    if cost==0: break        # Update feature matrix    hn=(transpose(w)*v)    hd=(transpose(w)*w*h)      h=matrix(array(h)*array(hn)/array(hd))    # Update weights matrix    wn=(v*transpose(h))    wd=(w*h*transpose(h))    w=matrix(array(w)*array(wn)/array(wd))        return w,h

参数pc指定我们希望找到的特征数，iter是最多执行次数。

来看看我机器执行的测试：

测试数据还是用之前的例子：每一行都是一本书，每一列都代表一个词，数组中的数字代表这个词在这本书中出现的次数或者TF-IDF值：

books = [      [2, 3, 3, 12, 13, 12],      [3, 3, 1, 13, 12, 11],   [1, 10, 3, 5, 11, 13],      [13, 12, 11, 1, 3, 3],      [12, 13, 12, 3, 2, 2],      [5, 11, 13, 3, 1, 10]  ]  

测试过程：

>>> from numpy import *>>> v=matrix(nnmf.books)>>> w,f=nnmf.factorize(v,2)2352.23524582130.308626905125.555272626125.540256481125.540046063>>> wmatrix([[ 1.01861466,  0.11646996],        [ 0.99224682,  0.07461746],        [ 0.80748193,  0.400123  ],        [ 0.04713232,  1.2421583 ],        [ 0.04119911,  1.27694878],        [ 0.27820173,  1.02855868]])>>> fmatrix([[  0.21984673,   3.29533888,   0.77713346,  10.7432405 ,          12.31672543,  12.3422249 ],        [  8.44995405,  10.06321664,   9.81092778,   0.6161926 ,           0.65151678,   2.88313112]])>>> w*fmatrix([[  1.20810496,   4.52874298,   1.93427795,  11.01499022,          12.62187925,  12.90776943],        [  0.84865631,   4.02068115,   1.50317469,  10.70592492,          12.26984613,  12.46166529],        [  3.5585432 ,   6.68745101,   4.55309905,   8.92152545,          10.20622012,  11.1197307 ],        [ 10.50654247,  12.65542508,  12.2233535 ,   1.27176259,           1.38980281,   4.16302294],        [ 10.79921601,  12.98597723,  12.56006946,   1.22945829,           1.33939163,   4.1900994 ],        [  8.75243536,  11.26737784,  10.30731484,   3.6225783 ,           4.09665751,   6.39909782]])

我测试用的特征数选择为2。

先来看看权值矩阵w：

book1，特征一明显(1.0)，特征二不明显(0.1)；

book2，特征一明显(0.9)，特征二不明显(0.0)；

book3，特征一明显(0.8)，特征二不明显(0.4)；

book4，特征一不明显(0.0)，特征二明显(1.2)；

book5，特征一不明显(0.0)，特征二明显(1.2)；

book6，特征一不明显(0.2)，特征二明显(1.0)；

再看看特征矩阵：

特征一：word1不明显(0.2)，word2不太明显(3.2)，word3不明显(0.7)，word4明显(10.7)，word5明显(12.3)，word6明显(12.3)

特征二：word1挺明显(8.4)，word2明显(10.0)，word3明显(9.8)，word4不明显(0.6)，word5不明显(0.6)，word6不明显(2.8)

实现还是很简单的，主要是这种思想，再次得成结论最NB的还是算法。上次听博士说一年要读100多篇论文，这才是学术研究，这才会创造核心技术。