白话NMF（Non-negative Matrix Factorization）——Matlab 实现
algorithmlixuan 2014-07-22
NMF——非负矩阵分解。如果你事先了解PMF[概率矩阵分解]的话，那么其实只要在PMF的基础上多加上一点，就是NMF了。
方法一：
在PMF中使用SGD【随机梯度下降】进行优化时，使用如下的迭代公式：
其中P、Q分别代表原始矩阵R的两个维度的隐含矩阵，在推荐应用中，一般讲P看做用户矩阵、Q看做物品矩阵。
从公式中不难看出，无论P矩阵还是Q矩阵都会出现负值的情况，上述公式并未对P、Q矩阵的值做任何限制。
在应用中，有时候需要分解出来的矩阵中不存在小于0的值，也即要求所有值非负。
怎么做到呢？其实很简单，在上述两个迭代公式中加个约束即可，如下公式：
很简单，是吧。这其实是非负矩阵分解实现中最常用的一种。
方法二：
在很早之前，大概是2001年的样子，Daniel D. Lee and H. Sebastian Seung.这两个家伙写了篇文章：
《Algorithms for non-negative matrix factorization》，其中讲了另外一种关于求解非负矩阵分解的方法，我们叫它迭代相乘法。
怎么做的呢？其实也很简单。
上一个方法中是用加减法来调整P、Q矩阵，既然加减不能保证非负，那用乘除是不是就可以？
如果我们将P、Q都初始化成非负矩阵，然后每次迭代都乘以一个非负的数（可以理解为“相对梯度”），这样是不是也可以呢？
当然是的。
具体可以这么理解：
- 如果某次的预测值比实际值大了，那么P，Q矩阵对应的值乘以一个小于小于1，并且大于0的数，
- 如果某次预测的值比实际值小了，那么P，Q矩阵对应的值就乘以一个大于1的数。
不过至于具体要乘的数是多大，则样看在迭代过程中，预测值和实际值之间的差距（相对梯度）来定。
迭代更新公式如下：
这个过程实现起来也极其方便，用matlab来帮助我们做一些矩阵乘法的体力劳动，代码如下：
%V为原始举证，W，H为分解之后得到的非负矩阵，R为降为之后的维数，K为迭代次数
function [W,H] = NMF(V,R,K)
[m,n]=size(V);
W = abs(rand(n,R));
H = abs(rand(R,m));
for i = 1:K
H = H .* (W'*V) ./ ((W'*W)*H);
W = W .* (V*H') ./ (W*(H*H'));
end
end假设我们进行如下测试：
m = rand(10,10); [w,h,l] = NMF(m,5,100);
结果为：
m：
0.9448 0.7722 0.9758 0.5674 0.4716 0.2537 0.8564 0.1323 0.5270 0.4794
0.7145 0.4754 0.5554 0.9688 0.5430 0.1326 0.8998 0.8705 0.8942 0.8985
0.6792 0.6809 0.8463 0.8245 0.0597 0.5450 0.2179 0.6030 0.7784 0.9347
0.9594 0.4169 0.4081 0.9596 0.6580 0.8278 0.0770 0.2653 0.0694 0.8179
0.7753 0.3801 0.4620 0.6463 0.8896 0.8370 0.4742 0.8648 0.2788 0.7089
0.6077 0.2133 0.8263 0.3796 0.1096 0.8333 0.8350 0.0581 0.3794 0.7432
0.9480 0.3829 0.9912 0.4766 0.4378 0.2037 0.4694 0.4578 0.8647 0.8997
0.0596 0.0297 0.5239 0.9119 0.2802 0.5444 0.4138 0.7222 0.4200 0.0652
0.2687 0.4723 0.9254 0.0149 0.9852 0.8749 0.5027 0.3390 0.2399 0.3359
0.9867 0.3334 0.7390 0.1567 0.6088 0.1210 0.1254 0.4012 0.5977 0.0043
w：
0.4435 0.9515 0.0007 0.1336 0.5970
0.2402 1.1113 1.2941 0.0026 0.0483
0.1836 0.9192 0.6268 0.8682 0.0030
0.7320 0.0007 0.4480 1.7279 0.0008
0.5884 0.0093 0.9411 0.9744 0.6506
0.0000 0.6284 0.0000 0.9964 0.6851
0.5200 1.0455 0.2000 0.3215 0.1344
0.0000 0.1706 1.3721 0.0005 0.5185
0.0969 0.0000 0.0156 0.4746 1.8222
1.0290 0.2060 0.0151 0.0000 0.2880
h：
0.8502 0.2643 0.3582 0.1457 0.5889 0.0000 0.0001 0.2576 0.2494 0.0295
0.4489 0.3566 0.5973 0.2882 0.0000 0.0003 0.5072 0.1090 0.5859 0.5611
0.0018 0.0040 0.0030 0.5250 0.1710 0.1694 0.1259 0.5310 0.1459 0.1062
0.1995 0.1082 0.0956 0.2491 0.0214 0.4638 0.0172 0.0000 0.0000 0.4622
0.0784 0.1731 0.4709 0.0000 0.4332 0.3818 0.3558 0.1370 0.0929 0.0025
w*h：
0.8776 0.5743 1.0210 0.3725 0.5228 0.2902 0.6974 0.3001 0.7237 0.6103
0.7097 0.4736 0.7767 1.0352 0.3837 0.2391 0.7439 0.8767 0.9043 0.7694
0.7433 0.4733 0.7011 0.8369 0.2352 0.5102 0.5611 0.4807 0.6761 0.9890
0.9682 0.3826 0.4295 0.7724 0.5451 0.8775 0.0868 0.4267 0.2484 0.8682
0.7515 0.3807 0.6187 0.8251 0.8102 0.8597 0.3715 0.7415 0.3499 0.5745
0.5346 0.4505 0.7933 0.4293 0.3182 0.7238 0.5796 0.1624 0.4319 0.8149
0.9864 0.5691 0.9054 0.5622 0.4056 0.2346 0.6089 0.3725 0.7840 0.7722
0.1198 0.1562 0.3502 0.7696 0.4593 0.4306 0.4438 0.8182 0.3483 0.2430
0.3199 0.3924 0.9381 0.1405 0.8593 0.9184 0.6584 0.2829 0.1957 0.2284
0.9899 0.3953 0.6272 0.2172 0.7333 0.1126 0.2089 0.3350 0.4063 0.1483
最后，我们看下这100次迭代中的w*h 与 m之间误差的变化情况：
方法三：
上面两中方法相对比较简单直接，还有一种方法是大名鼎鼎的 林智仁老爷子在2007年提出的。
《Projected Gradient Methods for Non-negative Matrix Factorization》，其中也对上述两种方法进行了介绍，不过主要是为了衬托自己方法的牛逼性。
其实该方法与第一种方法类似，只不过稍微复杂点：
- 首先，每次迭代对W、H不是指迭代一次，而是分别迭代多次，找到临时最优W、H。
- 其次，在每次子迭代时，通过一个孙子迭代寻找一个最优的步长alpha。
- 最后，林老爷子为NMF算法的退出，给出了一个理论的判断条件。
经过以上三点的改动，林老爷子在其文章中声称，他的方法比上述两种方式在收敛效率上有绝对的优势。
没错，你没听错，虽然步骤更多，过程更复杂，还在每步迭代中加了子迭代优化局部W、H，甚至在子迭代中还套了一个寻找最优步长alpha的孙子迭代。
但整体的收敛效率却更快了。
由于该方法比较复杂，感兴趣可以搜索林老的文章研读，同时林老爷子也给出了其方法的代码实现，这点要赞一下。
