非负矩阵分解NMF

非负矩阵分解NMF是一种新矩阵分解算法，自《Nature》1999年刊登了两位科学家D.D.Lee和H.S.Seung有关非负矩阵研究的成果后，此分解算法逐渐被人们接受并应用到各种领域。

    NMF的 基本思想可以简单描述为： 对于任意给定的一个 非负矩阵A，NMF算法能够寻找到 一个非负矩阵W和一个非负矩阵H，使得满足A=WH，从而将一个非负的矩阵分解为左右两个非负矩阵的乘积。

    NMF分解算法相较于传统的一些算法而言，具有实现上的简便性、分解形式和分解结果上的可解释性，以及占用存储空间少等诸多优点。通过矩阵分解，一方面将描述问题的矩阵的 维数进行削减，另一方面也可以对大量的数据进行 压缩和概括。在诸如图像处理、语义识别等都是 正值矩阵的应用，NMF提供了一种全新的分解方式，而且 分解结果往往带有实际意义。如对一个包含2429个脸部图像的数据库进行分解，可以得到结果正是人的鼻子、眼等人脸局部概念特征。事实上Lee和Seung在他们的论文中更深入地指出，与人类识别事物的过程相似，NMF也是一种优化的机制，近似于我们的脑分析和存储人脸数据的过程。

    更科普、更详细的介绍可以在网上搜《非负矩阵分解：数学的奇妙力量》，百度文库的链接如下 http://wenku.baidu.com/view/94c8af0bf78a6529647d5331.html 。

    值得注意的是，有人已经将 NMF应用于盲信号分离BSS领域，取得一定的效果，不过其研究的BSS有非欠定的要求（即源数不大于传感器数），且混合矩阵等都是正的。

    这里不过多介绍非负矩阵分解的理论，重点来说说它的应用和实现，实际上Matlab就有实现非负矩阵分解的函数 nnmf（至少2012b提供此函数）；除此之外，还有个 NMF信号/图像处理工具箱：NMFLAB for Signal and Image Processing，网址见 http://www.bsp.brain.riken.jp/ICALAB/nmflab.html ，不过貌似其只能在matlab2007b上运行，因为它只提供了2007b上编译的p文件； 163博客 http://fxy1211.blog.163.com/blog/static/68255322007826111015905/、 博客园博客 http://www.cnblogs.com/tiandsp/archive/2012/11/13/2768597.html有 matlab简单实现的代码，只有20行，能实现图像的压缩重构，如下图

    下面重点介绍一下matlab提供的nnmf函数及其用法，主要就是翻译一下matlab帮助了， 函数说明参见nnmf函数帮助， 用法举例可以在帮助里搜Feature Transformation，它里面介绍了几种特征变换，如NNMF、PCA、因子分解。

    nnmf函数说明：

[W,H] = nnmf(A,k)

    此函数是将nxm的非负矩阵A分解为非负矩阵W(nxk)和H(kxm)。分解的结果不是完全精确的， W*H只是A的一个低秩估计，这是因为分解算法的本质就是寻找W、H使得A与W*H之间的算术平方根残差D最小

D = sqrt(norm(A-W*H,'fro')/(N*M))

    分解的过程是 从随机初始W、H开始的迭代算法，由于残差优化目标D可能有局部最小，所以 不同的初始W、H可能会得到不同的结果。 如果算法的结果说收敛到一个比设定k更低的低秩，说明结果不是最优的，需要重新选择k值。

     函数得到的结果W、H做了一些标准化，如H的各行向量是单位长度的，W的各列向量是按长度降序排列的。

    函数的使用还可以指定额外的参数，如下

'algorithm'---

默认为'als'，表示使用交替最小二乘法；'mult'表示使用相乘更新算法。

通常'als'收敛的更快、更稳定，而'mult'对初始值更敏感，所以可以先用'mult'从多个随机初始值估计W、H，再用'als'从上面的W、H开始找更优结果。具体见下面应用。

'w0'---

设置W的初始值

'h0'---

设置H的初始值

'options'---

一个由statset函数创建的选项结构体，结构体内容包括：

Display — 显示的内容，'off'(默认) —不显示；'final'—显示最后结果；'iter'—显示迭代的各步结果。

MaxIter — 最大迭代数。默认是100，达到最大迭代数则当结果收敛对待。

TolFun — 迭代结束的残差容忍值。 Default is 1e-4.

TolX — 迭代结束的W、H容忍值。Default is 1e-4.

UseParallel — 'always'表示并行计算。Default is 'never'.

UseSubstreams — 'always'表示reproducible方式并行计算。Default is 'never'. 为了实现reproducibly式计算，需要设置Streams为下面的格式: 'mlfg6331_64' or 'mrg32k3a'.

Streams — 一个RandStream对象或其cell。如果不指定，将使用默认的stream产生随机初始值。如果指定了，除下面情况外使用单个对象：

打开了matlab pool或UseParallel是'always'或UseSubstreams是'never'。

'replicates'---

重复分解的次数，未指定w0和h0时全部使用随机初始值，否则也只是第一次使用设定w0、h0。它常和'mult'算法配合使用。默认是1。

    用法举例：

    先来看Feature Transformation里的例子吧，先来翻译matlab对非负矩阵分解的介绍：

    非负矩阵分解NMF是一个 基于特征空间低秩估计的降维技术，它在使特征维数降低的同时，保证特征都非负，从而蕴含一定的物理意义。

    分解的结果W的k列分别表示A中变量的变换，H的k行表示产生W中变换变量在A中n个原始变量中线性组合系数。k值的可能取值范围由实际情况确定。

    个人理解：上面的说法比较绕，可以从盲源信号的角度解释，如下图

    A就是X，W就是S。X是传感器测量得到的信号（各xi是列向量），S是实际的源，也就是说，那么多传感器测来测去都只是测那几个源的线性组合，所以传感器测得信号是冗余的，需要找到本来的几个源，于是就可以用NMF分解。从上面分解很容易理解， H的第一行是第一个源对各传感器的贡献，以此类推。所以各矩阵意义就很明确了。

    再来看一个例子，自带moore.mat文件中是生物化学中对氧气需求的5个预测器，加载信号并显示如下

load moore

X = moore(:,1:5);

 

figure

plot(X)

    下面使用 nnmf来对X做一个2秩估计，先使用'mult'算法从5个随机初始值来计算

opt = statset('MaxIter',10,'Display','final');

[W0,H0] = nnmf(X,2,'replicates',5,'options',opt,'algorithm','mult');

得到的结果如下

    rep       iteration       rms resid      |delta x|

      1          10         358.296      0.00190554

      2          10         78.3556     0.000351747

      3          10         230.962       0.0172839

      4          10         326.347      0.00739552

      5          10         361.547      0.00705539

Final root mean square residual = 78.3556

    可以看到5个初始值得到的结果都不一样，其中第4次的结果最好。下面 使用'als'方法从上面的结果出发做进一步优化

opt = statset('Maxiter',1000,'Display','final');

[W,H] = nnmf(X,2,'w0',W0,'h0',H0,'options',opt,'algorithm','als');

结果如下

    rep       iteration       rms resid      |delta x|

      1           2         77.5315     0.000830334

Final root mean square residual = 77.5315

    H的值如下

>> H

H =

    0.0835    0.0190    0.1782    0.0072    0.9802

    0.0559    0.0250    0.9969    0.0085    0.0497

     可以看到，X中第5个预测器（权重0.9802）对W中的第一个预测器影响最强；X中第3个预测器（权重0.9969）对W中的第2个预测器影响最强。

    为了更好的显示X中预测器的相对重要性，可以使用 biplot在W的列空间将结果画出来，如下

figure

biplot(H','scores',W,'varlabels',{'','','v3','','v5'});

axis([0 1.1 0 1.1])

xlabel('Column 1')

ylabel('Column 2')

图中也可以看到v5横坐标最大，即对W的第1预测器影响最大；v3的纵坐标最大，即对W的第2预测器影响最大。