非负矩阵分解(NMF)论文笔记（1）

NMF经典论文：Algorithms for Non-negative Matrix Factorization（NIPS, 2001），此文主讲算法细节偏理论。
（话说怎么去掉图片默认水印，有点影响视觉效果。）

1.引言

NMF要做的，就是给定一个非负矩阵$V\in\mathbb{R}_+^{d\times n}$，要找两个非负矩阵因子$W\in\mathbb{R}_+^{d\times c}$和$H\in\mathbb{R}_+^{c\times n}$，使得：

$$V\approx WH \tag{1}$$ $V$看作数据矩阵，其中每一列是一个$d$维的样例（sample）， $c$经常选取为比$d,n$都要小的一个数，事实上，NMF有很好的聚类特性， $c$可以看作是数据的聚类簇数。由于$W$和 $H$都比原始数据矩阵$V$要小，所以NMF也顺道实现了数据的压缩表示。

继续来看式(1)，现在我们仅仅关注一个样例，即$V$中的一列，用$v$表示，对应的$H$中的一列用$h$表示，则有：

$$v\approx Wh$$ 这个表述有什么意思呢？它表示呀， $W$的列可以看作是$d$维空间的 $c$个基，$h$指明了原始数据 $v$用这c个基可以如何表示—–$v$由W的列线性组合而成，组合系数就是 $h$中元素。$h$其实也就可以看作是原始数据 $v$的一种低维表达（d—>c）。因为所有的数据都要由这c个基来表达，所以要想式(1)或式(2)有更低的近似误差，这组基必须能很好的发现数据的隐含结构。

矩阵分解问题，在数值线性代数中，已经研究得很广泛，由于NMF是带有非负约束的矩阵分解，所以以前的一些方法并不能直接很好的运用。于是也有了这篇文章的工作。

2.损失函数

评估$V$和 $WH$间的近似程度，根据采用的不同误差度量准则，文章提出了两种损失函数。一个是平方欧式距离，一个是散度（divergence）。

如上图中最后一句话所述，当A,B矩阵满足元素和为1时，D(A||B)就是KL散度/相对熵。对应的，NMF的目标式就可以定义为如下两个：

$$\min_{W\geq 0, H\geq 0} \|V-WH\|_F^2 \\ \min_{W\geq 0, H\geq 0} D(V||WH)$$ 上面两个目标式，对W或者H中一者而言是凸的（即固定一个变量，优化另一个变量时是凸的），但是对W和H两者同时而言是非凸的。

3.乘法更新准则

求解上一节中带约束的优化问题，梯度下降法可能是最简单的方法，它容易实现，但是呢，收敛也比较慢。其它的方法，比如共轭梯度法收敛快，但是不易实现。此外，这些基于梯度的方法，都需要确定一个参数，即步长，而调整这个步长参数也是比较麻烦的。为此，在收敛速度和算法实现难易之间折中一下，就是文章中提出的乘法更新准则了。

这里都说要W和H在距离函数的驻点上（函数的一阶导为0是驻点，另，二阶导为0是拐点），应该是指函数对W和H要可导（因下一节中要求偏导，此处仍存疑）。

4.乘法更新准则Vs.加法更新准则

这一节是说，乘法更新准则其实可以由加法更新准则得到，算得上是加法更新准则的一种特殊情况，这里的加法更新准则，就是传统的梯度下降。下面，以平方欧式距离为例，我们用梯度下降来更新H。首先，求一下梯度：

$$\frac{\partial J}{\partial H} = -W^T(V-WH)$$ 所以，元素 $H_{ij}$可以这样更新（论文中的下标是 $a\mu$，这个无所谓）：
$$H_{ij} := H_{ij}-\eta_{ij} \frac{\partial J}{\partial H_{ij}} = H_{ij}+\eta_{ij}\left((W^TV)_{ij} - (W^TWH)_{ij}\right) \tag{6}$$ 这里步长 $\eta_{ij}$只要设置得足够小（大于0），式(6)的加法更新准则就可以保证目标函数值减小。

接下来，我们给步长取一个很特殊的值：

$$\eta_{ij}=\frac{H_{ij}}{(W^TWH)_{ij}} \tag{7}$$ 代入式(6)，即可得到上一节式(4)中关于H的更新准则。

关于散度准则的目标式，就直接贴图了。

但是这样给$\eta$取值，好像并不妥当，因为式(7)或者式(9)中的步长取值可能比较大，所以好像并不能保证由此得到的乘法更新准则一定能降低目标函数值。所以，在下一节，就是来证明这个问题的，我们将会看到，这样得到的乘法更新准则是能保证收敛的。

5.收敛性证明

收敛性证明是文章的一个重点，比较理论，但也正是最精妙处所在，证明过程将数学中的构造思想发挥得淋漓尽致。但是估计受nips篇幅限制，内容写的比较简洁，所以读起来也比较晦涩。下面边贴图，边补充论文中省略的推导过程。

我们还是以平方欧式距离为例进行说明，即下面对定理1进行证明，我们将会用到辅助函数，类似于EM算法中用到的，首先给出辅助函数的定义：

关于$G$，简单理解为，给定变量$h'$，$G(h,h’)$是关于变量$h$的函数。接下来是一个引理：

这个引理是说，给定$h^t$（迭代第t次时的变量值），我在最小化$G(h,h^t)$时，其实也实现了降低目标函数值F之目的。从下图，可以很清晰的看懂这个引理的证明过程。（proof中第一个小于等于是由于式(10)的左式，第三个等号是由于式(10)的右式，中间的小于等于是因为式(11).）

自自然然地，我们现在就想，NMF的目标式子，存不存在这样的辅助函数呢，如果存在，是什么？因为一旦能确定辅助函数，我们根据上面的定义和引理，就能证明前面定理的收敛性了。

文章很强势而直接地给出了辅助函数的表达式。贴图感受一下：

有时候，可以根据需要达到的目标，倒着推导我们达到目标需要哪些条件，我们可以怎样去构造，有时候，纯粹就是天才般的灵感，加上知识积累形成的直觉，直接就知道怎么做了。这里为什么要这么构造、能这么构造，关键还有这个K为什么这么取值，我目前能想到的只有，式(14)的形式是在跟目标式的二阶泰勒展开靠拢，而K跟步长的取值有关（倒数）。

下面，我们来证明引理2。在看论文中给出的证明过程前，我们先来对目标函数进行一些变换，注意我们将用$X_{i\cdot}$和$X_{\cdot j}$分别表示矩阵$X$的第i行和第j列，$X_{ij}$是第i行j列元素。

$$\begin{align} F(H) &= \frac{1}{2}\|V-WH\|_F^2 \\ &=\frac{1}{2}\sum_{j=1}^n\|V_{\cdot j}-WH_{\cdot j}\|_2^2 \quad (按列展开)\\ &=\frac{1}{2}\sum_{j=1}^n\sum_{i=1}^d(V_{ij}-W_{i\cdot}H_{\cdot j})^2 \\ &=\frac{1}{2}\sum_{j=1}^n\sum_{i=1}^d(V_{ij}-\sum_{a=1}^c W_{ia}H_{aj})^2 \quad (逐元素展开) \end{align}$$ 注意，这里F(H)是关于H的函数，在更新H时，变量W是固定不变的。上面推导的第二行是将矩阵范数按列展开写，这n个子式（对应n列）可以被解耦（是无关的），所以最小化F(H)，就等价于分别最小化这n个子式，我们以第j个子式为例写出新目标式：
$$F(H_{\cdot j})=\frac{1}{2}\sum_{i=1}^d(V_{ij}-\sum_{a=1}^c W_{ia}H_{aj})^2$$ 将V和H的第j列分别记为 $v$和$h$，则 $V_{ij}=v_i, H_{aj}=h_a$，上式可以简写为：
$$F(h)=\frac{1}{2}\sum_{i=1}^d(v_i-\sum_{a=1}^c W_{ia}h_a)^2$$ 至此，我们得到了引理2中的式(15)。我们对这个式子顺便求一下导吧：
$$\nabla F(h)=\frac{\partial F(h)}{\partial h} = -W^Tv+W^TWh\quad(对向量求一阶导，对按列展开的式子求比较简洁)$$ $$\frac{\partial F(h)}{\partial h\partial h} = W^TW\quad(对向量求二阶导，得到的即是Hessian矩阵)\\$$ 函数对向量/矩阵直接求导，是比较简洁明了的，如果不熟悉，也可以对各元素求偏导，然后把各元素求得的结果按顺序排放成向量/矩阵即可。
$$\frac{\partial F(h)}{\partial h_a} =\sum_{i=1}^d-W_{ia}(v_i-\sum_{a=1}^c W_{ia}h_a)=-W_{\cdot a}^Tv + W_{\cdot a}^TWh=-(W^Tv)_a+(W^TWh)_a\quad (对元素求一阶偏导，对式(15)求) \\ \frac{\partial F(h)}{\partial h_a\partial h_b} = \sum_{i=1}^d W_{ia}W_{ib}=W_{\cdot a}^TW_{\cdot b}\quad (对元素求二阶偏导，a,b\in[1,c])$$

接下来，我们就可以写出目标式F(h)的二阶泰勒展开：

再次把辅助函数表达式贴过来，对比着看：

很容易看出，G(h,h)=F(h)这个条件是满足的，下面我们只需要证明$G(h,h^t)\geq F(h)$，也即是要证明下式成立：

那么我们现在就是要证明矩阵$(K(h^t)-W^TW)$半正定，文章中用到了一个技巧，不是直接证它，而是考虑下面这个矩阵：

$M_{ab}$是对矩阵$(K(h^t)-W^TW)$中对应元素的一个缩放，故当且仅当M是半正定阵时，$(K(h^t)-W^TW)$也是。（关于这部分没有深究，还需要查资料多了解一下）

然后就是证明M半正定：

这里补充一下如何从式(19)到式(20)，把式(13)中的分子展开写：

$$(W^TWh^t)_a=(W^TW)_{a\cdot}h^t=\sum_b(W^TW)_{ab}h_b^t$$ 并注意有：
$$\sum_{a,b,a=b}=\sum_a$$
$$\begin{align} (19) &=\sum_{ab}\nu_ah_a^tK_{ab}(h^t)h_b^t\nu_b - \sum_{ab}\nu_ah_a^t(W^TW)_{ab}h_b^t\nu_b \\ &=\sum_{a,b,a=b}\nu_ah_a^t\frac{(W^TWh^t)_a}{h_a^t}h_a^t\nu_a - \sum_{ab}\nu_ah_a^t(W^TW)_{ab}h_b^t\nu_b \\ &=\sum_a \nu_a\sum_b(W^TW)_{ab}h_b^th_a^t\nu_a - \sum_{ab}\nu_ah_a^t(W^TW)_{ab}h_b^t\nu_b \\ &=(20) \end{align}$$

现在，我们已经证明了$G(h,h^t)$就是$F(h)$的辅助函数，所以我们可以通过最小化G()，来降低F()的值。令G对h的梯度（这里导数是向量，所以叫梯度）为0：

$$\frac{\partial G(h,h^t)}{\partial h} = \nabla F(h^t)+K(h^t)(h-h^t)=0$$ $$h^{t+1} \leftarrow h = h^t - K(h^t) ^{-1}\nabla F(h^t) \tag{24}$$
对每个元素，则有如下更新法则：
$$\begin{align} h_a^{t+1} &= h_a^t - K_{aa}(h^t)^{-1}\frac{\partial F(h)}{\partial h_a} \\ &= h_a^t - \frac{h_a^t}{(W^TWh^t)_a} \left(-(W^Tv)_a+(W^TWh^t)_a\right) \\ &=h_a^t\frac{(W^Tv)_a}{(W^TWh^t)_a} \tag{25} \end{align}$$
于是，我们推证了前面给出的关于H乘法更新准则，关于W的可以类似推证，此时目标式就应该是F(W)，具体过程以后有需要再推导。

关于散度目标函数情形，直接贴图了。