非负矩阵分解算法概述之Lee&Seung的世界.docx

非负矩阵分解算法概述
（吴有光）
NOTE：本文为科普文章，尽量做到通俗而不严格，比较适合理论小白补补NMF历史
第一部分 Lee&Seung的世界
1 引言
现实生活中的数据，我们总是希望有个稀疏表达，这是从压缩或数据存储的角度希望达到的效果。从另一方面来讲，我们面对大量数据的时候，总是幻想能够发现其中的“规律”，那么在表示或处理的时候，直接操作这些提纲挈领的“规律”，会有效得多。这个事情，让很多的科学家都伤透脑筋，不过也因此有了饭碗。
1.1 第一个例子
我们先来看一个简单的例子。在人文、管理或社会学里，实证研究方法是常用的方法。比如我们来考察大学生就业过程，对学生的选择工作类别的动机，我们常说“想吃劳保饭的同学铁了心要考公务员，喜欢轻松自由氛围的同学更趋向于外企，只想稳定的同学认为国企最好，富二代神马的最爱创业然后继承家产了”，这句话如果要严格来论证是不可能的，那么我们转而寻求“调查论证”，即通过设计问卷（问卷上设计了可能影响学生选择的因素，比如家庭情况、学业情况、性格取向、对大城市或家乡的热恋程度、以及人生观价值观等等各种我们可能会影响就业取向的因素）各种我们猜测会影响学生。
问卷上来后，我们通过统计得到如下的列表。

图1 第一个例子的统计表示例
表中的各个因素我们进行了量化，比如性格因素从完全内向到热情奔放分为5个等级（可以用一些问题来直接或间接获得这个等级）。那么剩下的问题就是回答开始的问题：
（1）是不是我们设计的每个因素都有效？（显然不是，之所以设计问卷就是要来解决这个问题的）
（2）是什么因素影响了学生的最终选择？或者说，从统计上来看，每个因素占多大比重？
这时，用矩阵来表示可写为，其中就表示那个因素矩阵，表示最终取向，代表我们要求的系数。我们把要求的用代替，写成矩阵形式为：
(1)
更进一步，如果我们不仅调查学生的去向，还想同时调查很多事情，那么就会有，这样上面的式子改写为：
(2)
此时问题转化为：
Q1：已知，如何求解，使之满足上面的等式，其中具有初始值（就是我们设计的一堆东西）。
如果我们让固定，这就是一个方程求解的过程。然而，当我们认为也可以缩减，即认为很少样本就足够表示我们真实取得的样本，那么问题进一步转化为：
Q2：如何同时求解和，使之满足。
或者我们也可以只对因素矩阵进行分解，即直接对其进行消减：
(3)
其中，为消减后因素矩阵，为在基底下的表示系数，这里要求列数要大大低于的列数，否则就没有实际意义。
上面这个过程，就类似Paatero&Tapper于1994年提出的实矩阵分解（Positive Matrix Factorization, PMF）模型，此模型后来被Lee&Seung提出的非负矩阵分解(Nonnegative Matrix Factorization, NMF/NNMF)模型所取代。
1.2 第二个例子
第一个例子为了给非数学、非信号处理的同学一个印象，写的罗里吧嗦，那第二个例子我们就简单写。
给定一组信号，如何找到对其进行稀疏表示？即如何找到满足的和，因为，这里要求且。
这个问题对信号处理的同学来说，太熟悉了。因为我们毕生的精力都在干这件事情。

如果去掉的非负限制，是有很多现成且高效的方法的，比如主成分分析（Principle Component Analysis， PCA）、独立成分分析（Independent Component Analysis， ICA）、因子分析（Factor Analysis， FA）等。然而，施加了非负限制后，这些方法就不适用了。而为什么要施加非负限制，回想第一个例子就明白了，我们最终找的是“影响因子”，因子会有负的么？

于是，非负矩阵分解就出世了，
1.3 非负矩阵分解
非负矩阵分解（Non-negative Matrix Factorization， NMF）从1999年正式提出【1】至今，经过十多年的发展，已经成为一个相对成熟的数据分析手段，在图像分析、文本聚类、数据挖掘、语音处理等方面得到了广泛应用。
NMF得到研究人员的青睐，除了易于获得快速分解算法之外，主要归功于其分解结果有较为明确的物理意义。例如在人脸识别中，分解结果为人脸的各个局部诸如鼻子、眼睛、嘴巴等，这符合人类思维中局部构成整体的概念。如图1所示

图1 Lee和Seung的经典文献中所使用的NMF说明图【1】
上图为NMF对人脸图像的分解结果，可见每一个子图都是人脸的某个局部；下左图为VQ分解结果，每一个子图就是某个原始样本；右下图为PCA分解结果，子图由特征脸和各级误差脸组成

据说，据Lee&Seung说，NMF由于在分解过程中做了非负限制，得到的结果会像图3上一样，每个子图（类似于基）是全图的一部分，这显然有别于我们往常所用的分解，并且更符合于人类直观视觉过程“局部组成整体”。
但是，我们的读者是勇于实践的，不管读不读完本文，他们都会自己coding或者到QQ群180291507共享里面下载代码来自己尝试尝试。
注释：在我试验的算法中，很难重现Lee&Seung的上述分解结果，在Naiyang Guan的最新NMF方法——MahNMF里【9】，也没有看到这种完全局部化的分解结果，不知道Lee&Seung是不是吹牛还是怎的。

下面从NMF问题描述出发，介绍NMF的发展历史和常见算法。
1.4历史渊源
矩阵分解是高维数据分析中重要方法，如主成分分析（Principle Component Analysis， PCA）、独立成分分析（Independent Component Analysis， ICA）、因子分析（Factor Analysis， FA）、矢量量化（Vector Quantization， VQ）等。而包含于各分析算法中的矩阵分解包括特征值分解、QR分解、LU分解等。关于这些分解，我在下一篇小白文里会介绍，等不及的筒子可以直接翻看矩阵论课本。
Paatero和Tapper于1994年尝试对环境数据进行因子分解时建立了如下优化模型【2】：
(4)
其中，为观测数据，为分解后左右矩阵，为权重矩阵。
三年后（1997年），Paatero将上述模型和求解方法进一步规范和推广为”Multilinear Engine”【3】，这时Paatero的算法已经现在乘性迭代的影子。但Paatero和Tapper的算法还存在一个大问题：算法收敛性和解的唯一性没有被证明！ 解唯不唯一倒是不重要，但算法的收敛性证明是重要的。
没证明啊没证明，没证明怎么行！
或许有人认为，好用就行了，管什么理论证明。如果这样想，就大错特错了。想想多么牛B的傅里叶变换是如何的历经磨难，由此导致傅里叶本人也是几经沉浮，唯一的原因就是傅里叶变换没有严格的理论证明，而当时的法国数学界是以严格出名的拉格朗日统治着。
现在和平时期，不搞战争不搞政变闲人多，我们立刻就能见到挥舞数学铲子跑过来填补漏洞的家伙。
这个数学铲子挥得好的第一人应该就是D. Lee和H. Seung（Seung现在是MIT脑认知科学实验室的头头），也是今天我们认为NMF的发明者【1】。他们把填漏洞的过程记录下来，顺手补了一个插图，扔给Nature杂志，于是就被录用了。于是他们这篇工程记录单就达到了三千多的引用率。这个故事告诉我们，你在家修灯泡修洗衣机的记录单或者心得体会都可以加几个数学公式，整巴整巴投到Science一不小心就中了。

图2 D.Lee 和 S. Seung（右）
我觉得他们俩是韩国人，元芳你怎么看？
后面还有很多挥铲子的，不过多是跟从者，包括因为写SVM代码包而大名鼎鼎的Chil-Jen Lin【4】，但在NMF界可就没人注意到他们的大名了。
1.2 问题描述
现在我们来看看Lee 和Seung对本问题的描述。从他们开始，这个问题就正式成为非负矩阵分解NMF。NMF是在所有矩阵元素非负的约束下的分解方式。形如：
(5)
式中为待分解矩阵，和为分解后两个矩阵。这里表示。
熟悉矩阵分解的同学可以看做是满秩分解。熟悉信号处理（正交分解，例如傅里叶变换）的同学可以把看做信号分解时候的基函数，而看做表示系数，在很多论文中作者也都是这么描述的，不过记住一点，基向量是相互正交的，而这里对矩阵的各向量无正交要求。如果换一种写法会更像以前学过的信号分解方式：
(6)
或者
(7)
式中为的第列，即第个基向量，为向量的第个元素，即表示系数。在傅里叶变换里会有
(8)
成立。
这里因为不是正交分解，所以不能使用式。
因为各向量不正交，那么就很难用构造的方式来获得（比如傅里叶变换矩阵，小波变化矩阵都是用构造得到的）。于是分解就面临如下问题：
（1） 已知如何同时确定？
（2） 分解中为数据降维的重要标志，如何选取？
针对问题（2），我们可以直接回答，的选取是人为的！ 也就是你自己一个一个试，比如，啦，你根据实际情况设置。不过可以想见的是原始数据为维，时，降维效果才好，数据量才能降下去。 但是，用脚趾头想问题都知道，太小了恢复效果肯定不好。
针对问题（1），就是本文的重点了，也是这么多人研究NMF的重点。
一个已知量，两个未知量怎么肯定求出来？？直接求，肯定没有办法，神仙也没办法。一些画瓢看多了的人就知道了，一定要加约束条件，加上约束条件问题就明朗了。
第一个约束条件是等式中隐含的。我们一般求解的问题中很难得到”=”，即我们得到的是近似解，于是式可以改写为：
(9)
其中是逼近误差。这样简单的一改，从理解问题的角度，可就大大的简单化了。
回头看式，要求我们找到两个矩阵，相乘等于，其他什么都不知道。这怎么找？再看式，要求我们找到两个矩阵，相乘近似等于，这一“近似”，可就有很多找法了。好比说，你找一个跟周杰伦一模一样的人来，你做不到；但要你找一个跟周杰伦有点像的人来，你就可以找一大堆。废话太多了，就是希望搞工程的筒子理解起来简单点，搞数值计算的人看家本领就是这样的近似来近似去的。

好了，到这里我们把问题重新描述一下，我们称为NMF问题。
NMF问题：已知，求解，使得
(10)
这里要求尽可能小，并且算法是快速收敛的。即找到求解下式的收敛算法
(11)
1.3 目标函数
在Lee 和Seung 的文章里，他们引入了两种目标函数，分别是欧式距离和K-L散度。以欧式距离来度量，式改写为
(12)
其中表示矩阵中的第行第列的那个元素。
如果以K-L散度为度量，则改写为
(13)
可以看出，只有当时，式和式才会得到最小值0.
似乎的似乎，问题并没有解决。如何同时求解式或式中的？做优化的告诉我们，对于或者，式是凸的，有最优解。但是对却不是，因为二者相乘了嘛。似乎这是个难题，不凸的问题不好求最优解。但是这种分解问题已经阻拦不了我们伟大的Lee和Seung了，他们从Paatero的文章思路出发，给出了交替迭代直至达到最优解的算法。
1.3 迭代规则
NOTE：下面这一段推导并不复杂，但只关心结果的同学可以直接看式、、和，这是求解和的迭代算法。
这里要感谢刘维湘同学，在看到Lee和Seung的原文之前，我是从他的论文中清晰的理解了NMF迭代公式的由来，推荐有兴趣的同学看一下文献【5】，保证你半小时看懂。
理解迭代规则之前，我们先回到式。要使最小，即求如下最大似然解【1,5,6】
(14)
一， 假设噪声服从高斯分布
则
(15)
式中。观察式和式，可将问题转化为求
(16)
的最小值。此时我们求出和即可迭代规则。
怎么迭代的？ 往负梯度方向走就是了。
这里要做一个假设，各数据点噪声的方差是一样的，即。这里不习惯矩阵乘法的同学先来熟悉一个式子：
(17)
所以
(18)
于是
(19)
(20)
其中为常数。
这样，我们就可从梯度负方向对实行迭代
(21)
(22)
再进一步，取，使两项抵消，得到乘性迭代规则:
(23)
(24)
二，假设噪声服从服从Poisson分布
这时，式改写为：
(25)
似然函数为
(26)
根据上面类似的推导，可以得到
(27)
(28)

至此，Lee&Seung的方法全部呈现完毕。我们总结为算法1：
算法1
已知数据举矩阵V和所能忍受的误差，求非负分解矩阵
（1） 随机初始化矩阵，要求非负；
（2） 应用迭代公式进行迭代，如果噪声为高斯噪声则根据式和式进行，如果噪声为Poisson噪声则根据式和式进行；
（3） 当误差小于时，或者达到最大迭代次数时，停止。

似乎，事情已经结束了，根据算法1去就可以达到我们的分解目的了。然而，有的同学会想到一个问题，算法1管用不？怎么来保证它管用？算法的更新是从负梯度方向进行的，那么如果能保证收敛性就好了，收敛性可以保证解的过程不会发生振荡。
好事做到底，Lee&Seung证明了算法1的收敛性【7】。我们下面给出简要过程。

1.4 收敛性证明
NOTE：算法收敛性证明是本文最为“数学”的部分，实在不感兴趣的同学请跳过。
艳阳高照的2001年，噢，或许是阴雨绵绵。Lee&Seung顺手把他们提出的两个迭代公式收敛性证明了【7】。
这个证明需要借助辅助函数，类似于EM算法（Expectation-Maximization，这是个好东西，他告诉我们了什么是交替迭代解非凸问题）中一样。既然是证明，自然少不了数学定义和数学推导，不喜欢的同学直接跳过此部分。
定义1 设是的辅助函数，如果下式满足
(29)
这么说有点抽样，我们来看个在教科书中被用烂了的图就明白了。

图3 辅助函数和关系图

引理1 如果是辅助函数，则在如下的更新序列下式非增的：
(30)
证明：，证毕。
引理2 如果是三角矩阵
(31)
则如下定义的是的辅助函数
(32)
(33)
证明： 显而易见。根据

将展开为
(34)
比较和的式子，可以将要证明的转化为

据说这是显然的【8】，但我没看出来。Lee和Seung使用半正定的定义对上式进行证明【7】，有兴趣的看看。
定理1 在算法1里面的式和式迭代序列下，是非增的。
证明： 对式中求并令为0得
(35)
其中为的极小值取到点。进而得
(36)
我们前面已经证明在此序列下式非增的，把和代入得
(37)
这个式子就和算法1里的迭代公式是一模一样的。

1.5 小结
第一部分重点介绍Lee和Seung的算法就到此结束了。如果你大概看懂了本文，并且按图索骥的看懂了文后所有参考文献的方法，那么第二部分就不用看了，因为它只是第一部分的延续而已，并无太多新意。
参考文献
【1】D. Lee, H. Seung; Learning the parts of objects by non-negative matrix factorization[J]. Nature, 1999, (1401):788-791
【2】P. Paatero, U. Tapper; Positive matrix factorization: A non-nagative factor model with optimal utilization of error estimates of data values[J]. Environ metrics, 1994,(15):111-126
【3】P. Paatero. Least-squares formulation of robust non-negative factor analysis[J]. Chemometrics and Intell. Lab. Sys. 1997, 1(37):23-35
【4】Chih-Jen Lin. Projected Gradient Methods for Non-negative Matrix Factorization[J]. 2007
【5】刘维湘, 郑南宁, 游屈波. 非负矩阵分解及其在模式识别中的应用[J]. 科学通报, 51(3):241-250, 2006
【6】P. Sajda, S. Du, L. Parra. Recovery of constituent spectra using non-negative matrix factorization[J]. Proc SPIE, 5207:321-331, 2003
【7】D. Lee, H. Seung. Algorithms for nonnegative matrix factorization[J], in Advances in Neural Information Processing 13, MIT Press, 2001, 556-562
【8】张宇飞. 加稀疏约束的非负矩阵分解[D]. 大连理工大学, 第12页
【9】Naiyang Guan, Dacheng Tao, Zhigang Luo, John Shawe-Taylo. MahNMF: Manhattan Non-negative Matrix Factorization[J], Journal of Machine Learning Research, submitted, 2012

有意交流者入群180291507，或围观博主
只是为了要代码的同学请不要加，因为一般不理睬上来就要代码的。NMF的代码网上也很多，也并不复杂，有兴趣了自己网上下载来看。
（未完待续）