CN104850533A - 一种带约束的非负矩阵分解方法及求解方法 - Google Patents

Abstract

本发明公开一种带约束的非负矩阵分解方法及求解方法，矩阵分解是将一个稀疏矩阵R分解为两个矩阵U和V的乘积的过程，从而用于预测稀疏矩阵R中未知值，因此广泛应用于各类推荐系统中。而非负矩阵分解在此基础上加上了U,V两个矩阵中所有元素均为非负的约束，使结果对于实际应用而言更具有可解释性。常见的非负矩阵分解方法往往较为复杂难以进一步增加其他约束条件，从而限制了其扩展性。本发明提出了一种易于扩展的，同时对矩阵R中未知值具有较高预测准确度的非负矩阵分解方法，并给出了其在推荐系统中应用的具体实施步骤。

Description

一种带约束的非负矩阵分解方法及求解方法

技术领域

本发明涉及矩阵分解算法，特别是侧重于非负矩阵分解算法，以及对应的基于随机梯度下降法的具体解法，并在此基础上提出了在推荐系统中的应用实例。

背景技术

矩阵在各类科学研究中都广为应用，矩阵分解可以将描述问题的原始矩阵的维数大大减少，起到压缩数据的作用，甚至还可以用来对原始数据中缺少的部分进行预测。

非负矩阵分解(NMF)是近年来提出的一种矩阵分解的新方法，已成功地应用于许多领域当中。较系统的NMF理论由D.D.Lee和H.S.Seung在1999年《Nature》上提出，并同时用于人脸识别，引起了科学界的广泛关注。

常见的传统矩阵分解方法有：主成分分析(PCA)、独立成分分析(ICA)、矢量量化(VQ)、奇异值分解(SVD)等。这些方法都允许结果中出现负值，这虽然在计算上是正确的，但是往往造成了结果的不可解释性，比如在下文中将提及的评分矩阵中，用户的评分就不可能是负数。而非负矩阵分解正是在矩阵分解过程中加上了所有元素均为非负这一限制条件，大大提高了结果的可解释性。

但是目前的NMF方法较为复杂，难以在其基础上进行扩展，如增加其他约束条件。

为此，我们提出了一种新的带约束的非负矩阵分解及求解方法。

发明内容

发明目的：由于现有的非负矩阵分解方法通常难以根据实际情况进行进一步的扩展，本发明提出了一种新的带约束的非负矩阵分解及求解方法，能在得到较高精度结果的同时，具有易于扩展的特效，能够有效地适用于各种实际的应用场景，并给出了其在推荐系统中的一种具体应用。

技术方案：一种带约束的非负矩阵分解方法及求解方法，包括如下部分：

1)带约束的非负矩阵分解算法；

2)基于随机梯度下降的分布式矩阵分解求解方法；

这种带约束的非负矩阵分解算法，区别于传统的NMF算法，由对R＝UV中 的U,V两个矩阵进行约束，来保证分解结果的非负性。设U’、V’为与U、V相同大小的两个矩阵，令U＝Sigmoid(U’)，V＝Sigmoid(V’)，其中Sigmoid(x)＝1/(1+e^-x)。由Sigmoid函数的定义可知，其定义域为(-∞,+∞)，而值域为(0,1)，由此可以保证U,V中的每一项都为非负值。但由于U,V中的项的取值范围被限定在(0,1)，因此还需要对R矩阵进行标准化，将其中的每个元素也转换到(0,1)上，才能保证分解的正确性。即令R_ij＝(R_ij-R_lower)/(R_upper-R_lower)，其中R_upper为R中元素取值的上界，R_lower为R中元素取值的下界。

在矩阵分解中，通常由损失函数来给出具体的解法。在本发明所提出的矩阵分解算法中，采用的是常见的平方损失函数：C＝(R-UV)²，此为最基础的形式。

在矩阵分解的具体求解过程中，通常会遇到过拟合的问题，即求得的模型对于用于训练的样本能得到很精确结果，但对于用来测试(或实际中需要预测)的样本却只能得到较差的结果。为了减轻过拟合的问题，通常在损失函数中加入正则项：C＝(R-UV)²+λ(‖U‖²+‖V‖²)。再代入对U,V的约束后可以得到C＝(R-S(U′)S(V′))²+λ(‖S(U′)‖²+‖S(V′)‖²)，其中S为Sigmoid函数。由公式argmin((R-S(U′)S(V′))²+λ(‖S(U′)‖²+‖S(V′)‖²))求得U’，V’，再由U＝Sigmoid(U’)，V＝Sigmoid(V’)最终求得U,V，即得到我们需要的模型。

公式argmin((R-S(U′)S(V′))²+λ(‖S(U′)‖²+‖S(V′)‖²))给出了U，V两个矩阵理论上的解法，本发明中采用随机梯度下降法来对该公式进行具体的求解。在实际求解中，R通常为一个稀疏矩阵，我们将这些已知的项称之为训练数据。要最小化损失函数C的值，实际上是要最小化C在训练集上的值，因此这里可以进一步展开为 $C={\mathrm{\Σ}}_{(i,j)\∈\mathrm{train}}{(R_{\mathrm{ij}}-{\mathrm{\Σ}}_{k=1}^{K}S\left(U_{\mathrm{ik}}\right)S\left(V_{\mathrm{jk}}\right))}^2+\mathrm{\λ}\left(\right||S\left(U_i\right){\left|\right|}^2+$ ${\left|\right|S\left(V_j\right)\left|\right|}^2,$ 其中i，j分别为评分中的用户序号和物品序号，λ为正则项系数，U_ik，V_jk，U_i就是表示矩阵里对应的项和行。在函数C中有两组参数(U_ik和V_jk)，对他们分别求偏导数可以得到：

$\∂C/{\∂U}_{\mathrm{ik}}=-2(R_{\mathrm{ij}}-U_i{V}_j)S\left(V_{\mathrm{jk}}\right)S^{\′}\left(V_{\mathrm{jk}}\right)+{2\mathrm{\λU}}_{\mathrm{ik}}$

$\∂C/{\∂V}_{\mathrm{jk}}=-2(R_{\mathrm{ij}}-U_i{V}_j)S\left(U_{\mathrm{ik}}\right)S^{\′}\left(U_{\mathrm{ik}}\right)+{2\mathrm{\λV}}_{\mathrm{jk}}$

S(x)＝1/(1+e^-x)

S′(x)＝e^x/(1+e^x)²

根据梯度下降法，需要将参数沿着最速下降方向向前推进，因此可以得到如下递推公式：

U_ik＝U_ik+α((R_ij-U_iV_j)S(V_jk)S′(V_jk)-λU_ik)

V_jk＝V_jk+α((R_ij-U_iV_j)S(U_ik)S′(U_ik)-λV_jk)

而随机梯度下降法，即是随机地抽取训练数据，用上述递推公式进行迭代，直到满足终止条件(在训练数据上的RMSE值的变动小于一个阈值)为止。

由以上描述可以看出，本发明所提出的矩阵分解方法的核心是其损失函数以及对应的最小化求解方法。而由于该损失函数形式上较为简单，可以根据不同应用场景的需要，较为便利地加入新的约束项，这就使得本方法具有了易于扩展的特性。

传统的梯度下降法并不支持并行计算，可扩展性差，难以应用于数据规模较大的场景。为此，本发明对所采用的随机梯度下降法进行了改良，使其可以支持分布式的计算环境。随机梯度下降法可以较为简单的修改为共享内存方式的并行算法，从而利用GPU之类的硬件进行快速计算，但受限于GPU内存容量，此类并行算法难以处理大规模的数据。将算法改造为分布式计算集群环境中的并行算法，则主要面临数据同步的开销问题，不同的节点在遇到训练数据时，可能会改变U、V矩阵中相同位置的元素值，这就造成了多个节点间数据冲突的问题，常见的解决方法是将训练数据分成N*N个块(N为节点数)，每一次在各个节点上分配互相独立的数据块进行计算，最后再进行一次U、V的同步，这也是现有的分布式随机梯度下降算法DSGD所采用的方法。而这样做所带来的问题是在一轮迭代中，需要进行N次同步。在分布式环境下，节点间的同步操作需要占用非常大的开销，因此我们希望能在此基础上减少同步次数，从而进一步提高算法的效率。

由于上述矩阵分解的优化问题不是一个凸优化问题，通常只能求得局部最优解。可以利用交替优化U、V矩阵的方法，寻找局部最优解。而将这种方法应用于我们的分布式随机梯度下降算法后，则有可能将每轮迭代的同步次数减少到1次。具体方法是，在奇数轮将数据按行分成N份(N为节点数)，在这一轮只对U矩阵进行更新，在这一轮完成之后，在所有节点间进行U矩阵的同步；在偶 数轮将数据按列分成N份，在这一轮只对V矩阵进行更新，在这一轮完成之后，在所有节点间进行V矩阵的同步。

在上述基础下，我们假设如下实际场景：某电影相关网站共记录有M部电影，N位用户，每位用户在网站上有若干条为电影进行评分的记录，现该网站需要一个为用户还未看过的电影的预测其评分。本发明使用矩阵分解的方法来实现这一功能。N位用户对M部电影的评分构成一个N*M的评分矩阵R，由于每位用户只分别对若干部电影有过评分，因此该R矩阵是一个稀疏矩阵。我们假设每一个用户均有自己对电影的偏好，用一个K维的向量u来表示，同时每一部电影有不同的特征，也用一个K维的向量v来表示，而这两个向量的内积就是用户对这部电影的评分。所有用户的偏好构成矩阵U，所有电影的特征构成向量V，写成矩阵乘法的形式即为R＝UV，恰好是上述矩阵分解方法中的分解形式。应用上述的矩阵分解方法，并用对应的随机梯度下降法进行求解，就可以得到U,V这两个矩阵，再利用等式UV＝R，可以得到完整的评分矩阵，其中原本就有的数据称为训练数据，原本空缺的数据称为预测值，预测值即所需要的数据。

有益效果：与现有技术相比，本发明提出一种兼具精度与高可扩展性的矩阵分解与求解方法。

附图说明

图1为评分数据的格式；

图2为评分数据示例；

图3为随机梯度下降法的算法流程图；

图4为分布式随机梯度下降法的算法流程图；

图5为一个实际推荐系统的实例，图中的模型训练即为本发明提出方法所应用的位置。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

1、输入输出数据说明

我们首先对应用中所涉及的输入输出数据，以及评价标准进行说明。

输入数据为用户对电影的评分矩阵，以稀疏矩阵的方式输入。每一行数据表示矩阵中的一项，标出了用户编号，电影编号，以及对应的评分。

输出数据为每个用户对所有电影的预测评分，标出了用户编号，电影编号及对应的预测评分。

2、初始化

在用随机梯度下降法进行矩阵分解的过程中，我们并不需要构建完整的R(M*N)矩阵，而只需要构造U(N*K),V(M*K)两个矩阵，这样可以节省大量的空间，K的取值一般可以在20～50之间。由于U,V矩阵用Sigmoid函数进行约束，所以实际上构造的是U’，V’两个矩阵。在初始化时，我们将这两个矩阵用0～1的随机数进行填充。

由于在本发明中所述的矩阵分解方法中，要求R矩阵中各项的取值范围在(0,1)上，因此要对输入数据进行标准化。以原始评分数据在1～5分为例：

r′＝(r-1)/4

其中r为原始评分，r’为标准化的评分。

3、模型训练 

为了加快收敛速度，我们对所采用的随机梯度下降法进行了修改，不再是随机地抽取训练数据，而是在每一轮迭代中以随机的顺序使用所有的训练数据。

对于训练集中的任一数据R_ij，对k从0到K以如下公式进行迭代：

U_ik＝U_ik+α((R_ij-U_iV_j)S(V_ik)S′(V_ik)-λU_ik)

V_jk＝V_jk+α((R_ij-U_iV_j)S(U_jk)S′(U_jk)-λV_jk)

在此为了进一步提高结果的精度，我们还对评分的预测值进行了平滑处理，公式为：

R＝F(UV)

F(x)＝1/(1+e^-7(x-0.5))

因此迭代公式进一步调整为：

U_ik＝U_ik+α((R_ij-F(U_iV_j))F′(U_iV_j)S(V_ik)S′(V_ik)-λU_ik)

V_jk＝V_jk+α((R_ij-F(U_iV_j))F′(U_iV_j)S(U_jk)S′(U_jk)-λV_jk)

F′(x)＝7e^7(x-0.5)/(1+e^7(x-0.5))²

如果是在分布式环境下，则可以使用分布式的梯度下降算法，在连续的两轮迭代中，依次按照相关公式只对U、V矩阵进行更新。

在分布式环境下，可以利用交替优化U、V矩阵的方法，寻找局部最优解。而将这种方法应用于我们的分布式随机梯度下降算法后，则有可能将每轮迭代的同步次数减少到1次。具体方法是，在奇数轮将数据按行分成N份(N为节点数)，在这一轮只对U矩阵进行更新，在这一轮完成之后，在所有节点间进行U矩阵的同步；在偶数轮将数据按列分成N份，在这一轮只对V矩阵进行更新，在这一轮完成之后，在所有节点间进行V矩阵的同步。

具体算法见图3，4。

4、输出结果 

在随机梯度下降法达到收敛后，我们得到了训练好的U,V矩阵，对于任意的评分r_ui，可以由公式计算得到，这里是两个向量的乘积,表示U中第i行乘以V中第j列。

表1是本发明方法与现有方法的对比结果，带星号的为本发明提出的方法。数据集为Movie-Len 100K，评价指标为均方根误差(RMSE)。

$\mathrm{RMSE}=\underset{r\∈S_{\mathrm{test}}}{\mathrm{\Σ}}\sqrt{\frac{{(r-r^{*})}^2}{\left|S_{\mathrm{test}}\right|}}$

其中S_test为测试数据集，r为评分的实际值，r^*为根据训练数据集分解得到的U、V矩阵的乘积所得到的评分的预测值。

如图5所示，一个实际推荐系统的示例，本文提出的方法属于其中模型训练的位置。

表1对比结果表

Method	RMSE
		GlobalAverage	1.117
UserAverage	1.036
		ItemAverage	0.983
SlopeOne	0.902
		UserItemBaseline	0.908
ItemKNNPearson	0.871
		FactorWiseMatrixFactorization	0.86
MatrixFactorization	0.857
		BiasedMatrixFactorization	0.854
＊Our Method	0.854

Claims (4)

1.一种带约束的非负矩阵分解方法及求解方法，其特征在于，包括如下部分：

1)带约束的非负矩阵分解算法；

2)基于随机梯度下降的分布式矩阵分解求解方法。

2.如权利要求1所述的带约束的非负矩阵分解方法及求解方法，其特征在于：

对R＝UV中的U,V两个矩阵进行约束，来保证分解结果的非负性；令U＝Sigmoid(U’)，V＝Sigmoid(V’)，其中Sigmoid(x)＝1/(1+e^-x)；由Sigmoid函数的定义可知，其定义域为(-∞,+∞)，而值域为(0,1)，由此可以保证U,V中的每一项都为非负值；但由于U,V中的项的取值范围被限定在(0,1)，因此还需要对R矩阵进行标准化，将其中的每个元素也转换到(0,1)上，才能保证分解的正确性。

3.如权利要求2所述的带约束的非负矩阵分解方法及求解方法，其特征在于：矩阵分解算法中，采用平方损失函数：C＝(R-UV)²；

在矩阵分解的具体求解过程中，在损失函数中加入正则项：C＝(R-UV)²+λ(‖U‖²+‖V‖²)；再代入对U,V的约束后可以得到C＝(R-S(U′)S(V′))²+λ(‖S(U′)‖²+‖S(V′)‖²)，其中S为Sigmoid函数；由公式argmin((R-S(U′)S(V′))²+λ(‖S(U′)‖²+‖S(V′)‖²))求得U’，V’，再由U＝Sigmoid(U’)，V＝Sigmoid(V’)最终求得U,V；

公式argmin((R-S(U′)S(V′))²+λ(‖S(U′)‖²+‖S(V′)‖²))给出了U，V两个矩阵理论上的解法，采用随机梯度下降法来对该公式进行具体的求解；在实际求解中，R通常为一个稀疏矩阵，我们将这些已知的项称之为训练数据；要最小化损失函数C的值，实际上是要最小化C在训练集上的值，因此这里可以进一步展开为 $C={\mathrm{\Σ}}_{(i,j)\∈\mathrm{train}}{(r_{\mathrm{ij}}-{\mathrm{\Σ}}_{f=1}^{F}S\left(U_{\mathrm{ik}}\right)S\left(V_{\mathrm{jk}}\right))}^2+\mathrm{\λ}({\left|\right|S\left(U_i\right)\left|\right|}^2+{\left|\right|S\left(V_j\right)\left|\right|}^2;$ 在函数C中有两组参数(U_ik和V_jk)，对他们分别求偏导数可以得到：

$\∂C/{\∂U}_{\mathrm{ik}}=-2(R_{\mathrm{ij}}-U_i{V}_j)S\left(V_{\mathrm{jk}}\right)S^{\′}\left(V_{\mathrm{jk}}\right)+2\mathrm{\λ}{U}_{\mathrm{ik}}$

$\∂C/{\∂V}_{\mathrm{jk}}=-2(R_{\mathrm{ij}}-U_i{V}_j)S\left(U_{\mathrm{ik}}\right)S^{\′}\left(U_{\mathrm{ik}}\right)+2\mathrm{\λ}{V}_{\mathrm{jk}}$

S(x)＝1/(1+e^-x)

S′(x)＝e^x/(1+e^x)²

根据梯度下降法，需要将参数沿着最速下降方向向前推进，因此可以得到如下递推公式：

p_ik＝p_ik+α((R_ij-U_iV_j)S(V_jk)S′(V_jk)-λU_ik)

q_jk＝q_jk+α((R_ij-U_iV_j)S(U_ik)S′(U_ik)-λV_jk)

而随机梯度下降法，即是随机地抽取训练数据，用上述递推公式进行迭代，直到满足终止条件为止。

4.如权利要求2所述的带约束的非负矩阵分解方法及求解方法，其特征在于：对所采用的随机梯度下降法进行了改良，使其可以支持分布式的计算环境；在分布式环境下，利用交替优化U、V矩阵的方法，寻找局部最优解：在奇数轮将数据按行分成N份(N为节点数)，在这一轮只对U矩阵进行更新，在这一轮完成之后，在所有节点间进行U矩阵的同步；在偶数轮将数据按列分成N份，在这一轮只对V矩阵进行更新，在这一轮完成之后，在所有节点间进行V矩阵的同步。