非负矩阵分解(NMF)是1999年由Lee和Seung^[1]首次提出的, 有着比较广泛的应用, 如图像处理^[1]、文本聚类^[2]和社区发现^[3]等.相较于其他的矩阵分解, 由于元素的非负性, 所以可解释性很强.比如:在处理图像时, 负的像素点是没有意义的; 在文本聚类时, 正值代表文档以一定概率属于某个主题, 负值则没有任何实际意义.

非负矩阵分解数学描述如下:给定非负矩阵$ A\in \mathbb{R}^{m\times n} $, 希望找到两个非负矩阵$ W\in \mathbb{R}^{m\times k} $和$ H\in \mathbb{R}^{k\times n} $, 使得$ A\approx WH $.可以通过以下优化问题来找到$ W $和$ H $,

$ \min \limits_{W, H\geq 0} f(W, H) = \frac{1}{2}\| A-WH \| _F^2, $

其中$ \| \cdot \|_F $代表矩阵的Frobenius范数.

解决非负矩阵分解的方法有很多, 比较常见的是由Lee和Seung在2001年提出的基准算法(Multiplicative Update Rules)^[4]、Lin在2007年提出的投影梯度法^[5]以及Kim和Park在2008年提出的交替非负最小二乘法(ANLS)^[6].

不同于某些迭代算法, 非负矩阵分解初始值的选择对于算法效果有很大影响.由Qiao提出的SVD-NMF^[7]和Boutsidis等人提出的NNDSVD^[8]算法都是基于矩阵$ A $的奇异值分解, 但是当矩阵$ A $的维数很大时, 对$ A $直接进行奇异值分解是耗时的.不同于这两种方法, 本文不直接对原矩阵$ A $进行奇异值分解.本文利用Frieze等人^[9]的思想, 通过抽样构造一个更小的矩阵, 对这个小矩阵来做奇异值分解, 从而实现对$ W $和$ H $的初始化, 这大大节省了算法的计算时间.

符号说明:用$ \| A\|_{F} $表示矩阵$ A $的Frobenius范数; $ \| x \| $表示向量的2 -范数; 矩阵$ A $的第$ i $行和第$ j $列分别用列向量$ A(i, \cdot) $和$ A(\cdot, j) $来表示; $ A_{ij} $和$ x(i) $分别表示矩阵$ A $的第$ (i, j) $个元素和向量$ x $的第$ i $个元素; $ A^{\mathrm{T}} $和$ x^{\mathrm{T}} $分别表示矩阵$ A $和向量$ x $的转置.对于一个非零向量$ x\in \mathbb{R}^{n} $, 用$ D_{x} $表示$ x $在$ \{1, 2, \ldots, n\} $上的分布, 其中概率密度函数为$ D_{x}(i) = |x(i)|^{2}/\| x\|^{2} $.通常把依据分布$ D_{x} $产生的样本称为从$ x $中产生的样本.

首先, 介绍一种最为经典的非负矩阵分解算法, 即由Lee和Seung在2001年提出的基准算法(Multiplicative Update Rules)^[4], 以下简称MU算法.在本文的数值实验部分, 将会用到MU算法. MU算法步骤如下

步骤1 初始化$ W_{ia}^0 \geq 0 $, $ H_{bj}^0 \geq 0 $, $ \forall i, a, b, j. $

步骤2 对于$ k = 1, 2, \cdots $,

$ \begin{align} & W_{ia}^{k} = W_{ia}^{k-1} \frac{(A(H^{k-1})^{\mathrm{T}})_{ia}}{(W^{k-1}H^{k-1}(H^{k-1})^{\mathrm{T}})_{ia}}, \quad \forall i, a, \end{align} $

(2.1)

$ \begin{align} & H_{bj}^{k} = H_{bj}^{k-1} \frac{((W^{k})^{\mathrm{T}}A)_{bj}}{((W^{k})^{\mathrm{T}}W^{k}H^{k-1})_{bj}}, \quad \forall b, j. \end{align} $

(2.2)

从式子(2.1)和(2.2)中可以看出, 如果在第$ k $次迭代中$ W $和$ H $是非负的, 那么在第$ k+1 $次迭代中, MU算法依然能够保证$ W $和$ H $是非负的.并且, Lee和Seung^[4]已经证明在MU算法下, 目标函数$ f(W, H) $是非增的, 即$ f(W^{k}, H^{k-1})\leq f(W^{k-1}, H^{k-1}) $, $ f(W^{k}, H^{k})\leq f(W^{k}, H^{k-1}) $.

由于本文提出的算法涉及到抽样, 下面给出对向量和矩阵进行抽样的假设.

假设2.1^[10] 给定向量$ x\in \mathbb{R}^{n} $, 假设可以根据分布$ D_{x} $来对$ i\in\{1, 2, \ldots, n\} $进行抽样, 也就是说$ i $被抽中的概率为$ D_{x}(i) = \frac{| x(i)|^{2}}{\| x\|^{2}}. $

假设2.2^[10] 给定矩阵$ A\in \mathbb{R}^{m\times n} $, 假设满足以下假设

1.可以对矩阵$ A $的行标$ i\in\{1, 2, \cdots, m\} $进行取样, 矩阵$ A $的第$ i $行被选中的概率为

$ P_{i} = \frac{\| A(i, \cdot)\|^{2}}{\| A\|_{F}^{2}}. $

2.对于任意的$ i\in\{1, 2, \cdots, m\} $, 可以根据分布$ D_{A(i, \cdot)} $来对矩阵$ A $的列标$ j\in\{1, 2, \cdots, n\} $进行抽样, 也就是说$ j $被抽中的概率为

$ D_{A(i, \cdot)}(j) = \frac{| A_{ij}|^{2}}{\| A(i, \cdot)\|^{2}}. $

定理2.3^[9] 给定矩阵$ A\in \mathbb{R}^{m\times n} $, 独立的按照概率分布

$ \{\| A(1, \cdot)\|^{2}/\| A\|_{F}^{2}, \cdots, \| A(m, \cdot)\|^{2}/\| A\|_{F}^{2}\}, $

从行标中取$ p $个样本: $ i_{1}, \cdots, i_{p} $.记$ S\in \mathbb{R}^{p\times n} $为规范化后的子矩阵, 即对于$ t\in\{1, \cdots, p\} $, $ S(t, \cdot) = A(i_{t}, \cdot)/\sqrt{p\| A(i_{t}, \cdot)\|^{2}/\| A\|_{F}^{2}}. $则对于任意给定的$ \theta>0 $, 满足

$ {\rm Pr}(\| A^{\mathrm{T}}A-S^{\mathrm{T}}S\|_{F}\geq\theta\| A \| _{F}^{2})\leq \frac{1}{\theta^{2}p}. $

在第3节, 会根据定理2.3说明本文提出的初始化方法(KFV-NMF)的可行性.

由定理2.3, 可以看出, 若按照满足假设2.2的概率分布来抽样:假设抽取了矩阵$ A $的$ p $行, 并按照定理2.3所说的构造一个$ p \times n $阶的矩阵$ S $, 则可以用$ S^{\mathrm{T}}S $来近似$ A^{\mathrm{T}}A $. $ A $的奇异值分解可以通过$ A^{\mathrm{T}}A $的谱分解得到, 尽管这里可以从$ S^{\mathrm{T}}S $的谱分解入手, 但并不打算这样做, 因为$ S^{\mathrm{T}}S $是一个$ n \times n $阶的矩阵, 维数也是很大的.这里可以对$ S $的列进行同样的操作, 从$ S $中抽取$ p $列, 类似的得到一个$ p \times p $阶的矩阵$ W $, 则$ WW^{\mathrm{T}} $近似$ SS^{\mathrm{T}} $, 故最终可以直接对一个$ p \times p $阶的矩阵$ W $做奇异值分解, 从而近似得到矩阵$ A $的右奇异向量, 这大大减少了计算时间.这种思想是Frieze等人^[9]提出的, 本文借鉴这种思想提出了一种更节时的非负矩阵分解初始化方法, 记此算法为KFV-NMF.

下面, 给出算法的具体流程.

下面, 给出一个定理, 并以此来说明算法FKV-NMF中第2步的具体实施方法.

定理3.1 假设矩阵$ A\in \mathbb{R}^{m\times n} $满足假设2.2, $ S $和$ W $如算法FKV-NMF中第1和2步所构造, 则$ \| A \| _F^2 = \| S \| _F^2 = \| W \| _F ^2. $

证 由$ P_{i_t} = \| A(i_t, \cdot)\|^{2}/\| A\|_{F}^{2} $, $ S(t, \cdot) = A(i_{t}, \cdot)/\sqrt{pP_{i_t}} $, 有

$ \| S \| _F^2 = \sum\limits_{t = 1}^p \| S(t, \cdot) \| ^2 = \sum\limits_{t = 1}^p \frac{\| A(i_t, \cdot) \| ^2}{pP_{i_t}}\\ = \sum\limits_{t = 1}^p \frac{\| A \| _F^2}{p} = \| A \| _F^2. $

接下来证明$ \| S \| _F^2 = \| W \| _F ^2 $.

$ \begin{align*} P_j'& = \sum\limits_{t = 1}^p\frac{D_{A(i_t, \cdot)}(j)}{p} = \sum\limits_{t = 1}^p \frac{| A(i_t, j)| ^2}{p \| A(i_t, \cdot) \|^2} = \frac{1}{\| A \|_F^2}\sum\limits_{t = 1}^p \frac{\| A \| _F^2}{\| A(i_t, \cdot) \| ^2}\frac{| A(i_t, j)| ^2}{p}\\ & = \frac{1}{\| A \|_F^2}\sum\limits_{t = 1}^p | S(t, j) | ^2 = \frac{\| S(\cdot, j) \| ^2}{\| A \| _F^2} = \frac{\| S(\cdot, j)\|^2}{\| S \|_F^2}. \end{align*} $

由$ P_{j_t}' = \| S(\cdot, j_t)\|^{2}/\| S\|_{F}^{2} $, $ W(\cdot, t) = S(\cdot, j_{t})/\sqrt{pP_{j_t}'} $, 有

$ \| W \| _F^2 = \sum\limits_{t = 1}^p \| W(\cdot, t) \| ^2 = \sum\limits_{t = 1}^p \frac{\| S(\cdot, j_{t}) \| ^2}{pP_{j_t}'}\\ = \sum\limits_{t = 1}^p \frac{\| S \parallel _F^2}{p} = \| S \| _F^2. $

由定理3.1可知, 在算法FKV-NMF的第2步中定义的$ P_{j}' $, 经过计算, 即为$ \| S(\cdot, j)\|^{2}/\| S\|_{F}^{2} $, 故$ W $的取法和$ S $的取法本质上是一样的, 这也为数值实验带来了便利.但这并不说明算法第2步中$ P_{j}' $的取法是无意义的, 因为当矩阵$ A $很大时, 为了节省内存, 在知道抽样后得到的行标$ i_1, i_2, \ldots, i_p $后, 可以通过行标来访问原矩阵$ A $中元素, 从而间接得到$ S $的信息, 而不是把$ S $显式计算出来.

由算法FKV-NMF的流程可以看出, 本文提出的初始化算法FKV-NMF的数值计算时间大大减少的原因在于:本文不是对原矩阵$ A $直接进行奇异值分解, 而是通过抽样构造一个更小的矩阵$ W $, 对小矩阵$ W $做奇异值分解, 这使得算法FKV-NMF在奇异值分解上花费的时间大大减少, 而现有的基于奇异值分解的初始化算法都是直接对原矩阵$ A $进行奇异值分解, 这是本文提出的算法数值计算时间较其他算法减少的原因.另外, 通过简单计算可以发现:算法FKV-NMF中矩阵$ \hat{V} $中的每一列可以近似看作矩阵$ A $的右奇异向量, Frieze等人在其文章^[9]中也说明了这一点, 这也是算法FKV-NMF的理论依据.

当然, 算法FKV-NMF也有弊端:此算法虽然在奇异值分解上节省了很多时间, 但是它还要计算用于抽样的概率分布$ P $和$ P' $, 这是有点耗时的.一般来说, 基于奇异值分解的初始化算法使用的都是截断奇异值分解, 比如前文提到的初始化算法SVD-NMF和NNDSVD.随着$ k $值的减小, 截断奇异值分解所使用的时间也在不断减少, 当$ k $小到某一程度时, 算法FKV-NMF在进行奇异值分解上节省的时间便不再占优势.但是, 这也是算法FKV-NMF可以进一步改进的地方, 比如可以考虑使用其他更加快速有效的抽样方法.

在本次实验中, 测试了两组数据:第一组是随机生成的服从均值为$ 0 $、方差为$ 1 $的高斯分布的随机数, 这里取绝对值来保证矩阵元素的非负性, 即$ A = | N(0, 1) | $, 维数是$ 500 \times 300 $.第二组是人脸数据库ORL, ORL数据库一共包含$ 400 $张$ 92 \times 112 $像素的图片, 一共$ 40 $个人, 每人$ 10 $幅不同表情的脸.在计算不同算法的时间和误差的实验中, 将这$ 400 $张人脸图像组成了一个$ 10304 \times 400 $的矩阵, 作为算法的输入.在进行人脸重构实验时, 选取其中$ 360 $张图像作为训练集, $ 40 $张图像作为测试集.

在实验中, 将本文提出的算法FKV-NMF与两种比较常见的基于SVD的初始化算法SVD-NMF和NNDSVD作了比较.非负矩阵分解算法采用的是经典的MU算法.误差采用公式$ {\rm error} = \| A-WH \|_F /\| A \| _F $.

在表 1和表 2中, 分别计算了在$ k $等于不同值时, SVD-NMF、NNDSVD和FKV-NMF三种不同的初始化算法所用的时间和所产生的误差.另外, 由于本文提出的算法涉及到抽样技术, 所以算法效果会有稍许波动, 为公平起见, 实验重复20次, 取平均值.

从表 1和表 2中, 可以看出:就时间来方面说, 本文提出的FKV-NMF初始化算法所用的时间要比SVD-NMF和NNDSVD小的多.就误差方面来说, 初始化误差比NNDSVD要大一些, 但是比SVD-NMF要小.所以本文提出的初始化算法在大大节约了时间的同时, 误差也在可接受的范围内.

在图 1中, 展现了SVD-NMF、NNDSVD和FKV-NMF三种不同初始化算法对于非负矩阵分解算法误差的影响.图中曲线表示随着迭代次数增加, 误差的变化情况.从图 1中可以看出, 在使用SVD-NMF、NNDSVD和FKV-NMF三种不同的初始化算法后, 使用MU算法进行非负矩阵分解时, 随着迭代次数的增加, 误差均在减少, 直到达到一个平稳值. NNDSVD的初始化误差虽然最小, 但误差很快便不再下降.在随机生成的数据(Random)上, 大约100次迭代后, 误差几乎不再下降, 在ORL数据集上, 大约200次迭代后, 误差几乎也不再下降, 而且最终误差比其他两种方法都要高的多.本文提出的FKV-NMF虽然最终误差比SVD-NMF稍高, 但是相差并不多.而且在迭代的最初始阶段, 使用本文提出的FKV-NMF初始化方法的误差要比SVD-NMF低一些.

下面, 展示了使用不同初始化算法后, MU算法迭代次数为$ 1000 $时, ORL数据库部分人脸图像的重构情况. 图 2是人脸的重构图像, 第1到第3行使用的初始化算法依次是SVD-NMF、NNDSVD和FKV-NMF.在表 3中, 计算了重构图像的信噪比(SNR), 第1到第7列分别表示7个重构图像.信噪比是用来衡量图像质量的一种方式, 信噪比越大, 表示重构的人脸图像质量越好.这里信噪比的计算公式为: $ {\rm SNR} = 10*{\rm log}_{10}(\sum\limits_{i, j}g(i, j)^2/\sum\limits_{i, j}(f(i, j)-g(i, j))^2) $, 其中$ g(i, j) $表示原始图像的灰度值, $ f(i, j) $表示重构图像的灰度值.

从图 2和表 3中, 可以看出:从重构图像的清晰度来看, 使用SVD-NMF和本文提出的FKV-NMF初始化算法图像清晰度相差无几, 而使用NNDSVD初始化算法的重构图像效果则稍差.从SNR的值来看, 使用SVD-NMF和本文提出的FKV-NMF初始化算法图像的SNR值相差并不多, 而使用NNDSVD初始化算法的重构图像SNR值要明显小于其他两者.

通过以上数值实验, 可以看到:本文提出的初始化算法FKV-NMF的运行时间较SVD-NMF和NNDSVD大大减少, 也在一定程度上保持了精度.