推荐系统算法介绍

SVD:

 

SVD算法详解

下面开始介绍SVD算法,假设存在以下user和item的数据矩阵:

这是一个极其稀疏的矩阵,这里把这个评分矩阵记为R,其中的元素表示user对item的打分,“?”表示未知的,也就是要你去预测的,现在问题来了:如何去预测未知的评分值呢?上一篇文章用SVD证明了对任意一个矩阵A,都有它的满秩分解:

那么刚才的评分矩阵R也存在这样一个分解,所以可以用两个矩阵P和Q的乘积来表示评分矩阵R:

上图中的U表示用户数,I表示商品数。然后就是利用R中的已知评分训练P和Q使得P和Q相乘的结果最好地拟合已知的评分,那么未知的评分也就可以用P的某一行乘上Q的某一列得到了:

这是预测用户u对商品i的评分,它等于P矩阵的第u行乘上Q矩阵的第i列。这个是最基本的SVD算法,那么如何通过已知评分训练得到P和Q的具体数值呢?

假设已知的评分为:

则真实值与预测值的误差为:

继而可以计算出总的误差平方和:

只要通过训练把SSE降到最小那么P、Q就能最好地拟合R了。那又如何使SSE降到最小呢?下面介绍一个常用的局部优化算法。

梯度下降法

为了说明梯度下降法,我找了一张PPT:

也就是说如果要最小化目标函数,必须往其负梯度方向搜索。这就是梯度下降法,注意它是一个局部优化算法,也就是说有可能落到局部最优解而不是全局最优解。

Basic SVD

利用梯度下降法可以求得SSE在Puk变量(也就是P矩阵的第u行第k列的值)处的梯度:

利用求导链式法则,e^2先对e求导再乘以e对Puk的求导:

由于

所以

上式中括号里的那一坨式子如果展开来看的话,其与Puk有关的项只有PukQki,其他的无关项对Puk的求导均等于0

所以求导结果为:

所以

为了让式子更简洁,令

这样做对结果没有影响,只是为了把求导结果前的2去掉,更好看点。得到

现在得到了目标函数在Puk处的梯度了,那么按照梯度下降法,将Puk往负梯度方向变化:

令更新的步长(也就是学习速率)为

则Puk的更新式为

同样的方式可得到Qik的更新式为

得到了更新的式子,现在开始来讨论这个更新要怎么进行。有两种选择:1、计算完所有已知评分的预测误差后再对P、Q进行更新。2、每计算完一个eui后立即对Pu和qi进行更新。这两种方式都有名称,分别叫:1、批梯度下降。2、随机梯度下降。两者的区别就是批梯度下降在下一轮迭代才能使用本次迭代的更新值,随机梯度下降本次迭代中当前样本使用的值可能就是上一个样本更新的值。由于随机性可以带来很多好处,比如有利于避免局部最优解,所以现在大多倾向于使用随机梯度下降进行更新。

RSVD

上面就是基本的SVD算法,但是,问题来了,上面的训练是针对已知评分数据的,过分地拟合这部分数据有可能导致模型的测试效果很差,在测试集上面表现很糟糕。这就是过拟合问题,关于过拟合与欠拟合可以看一下这张图

第一个是欠拟合,第二个刚好,第三个过拟合。那么如何避免过拟合呢?那就是在目标函数中加入正则化参数(加入惩罚项),对于目标函数来说,P矩阵和Q矩阵中的所有值都是变量,这些变量在不知道哪个变量会带来过拟合的情况下,对所有变量都进行惩罚:

这时候目标函数对Puk的导数就发生变化了,现在就来求加入惩罚项后的导数。

括号里第一项对Puk的求导前面已经求过了,第二项对Puk的求导很容易求得,第三项与Puk无关,导数为0,所以

同理可得SSE对qik的导数为

将这两个变量往负梯度方向变化,则更新式为

这就是正则化后的SVD,也叫RSVD。

加入偏置的SVD、RSVD

关于SVD算法的变种太多了,叫法也不统一,在预测式子上加点参数又会出来一个名称。由于用户对商品的打分不仅取决于用户和商品间的某种关系,还取决于用户和商品独有的性质,Koren将SVD的预测公式改成这样

第一项为总的平均分,bu为用户u的属性值,bi为商品i的属性值,加入的这两个变量在SSE式子中同样需要惩罚,那么SSE就变成了下面这样:

由上式可以看出SSE对Puk和qik的导数都没有变化,但此时多了bu和bi变量,同样要求出其更新式。首先求SSE对bu的导数,只有第一项和第四项和bu有关,第一项对bu的求导和之前的求导类似,用链式法则即可求得,第四项直接求导即可,最后可得偏导数为

同理可得对bi的导数为

所以往其负梯度方向变化得到其更新式为

这就是修改后的SVD(RSVD)。

ASVD

全称叫Asymmetric-SVD,即非对称SVD,其预测式子为

R(u)表示用户u评过分的商品集合,N(u)表示用户u浏览过但没有评过分的商品集合,Xj和Yj是商品的属性。这个模型很有意思,看预测式子,用户矩阵P已经被去掉了,取而代之的是利用用户评过分的商品和用户浏览过尚未评分的商品属性来表示用户属性,这有一定的合理性,因为用户的行为记录本身就能反应用户的喜好。而且,这个模型可以带来一个很大的好处,一个商场或者网站的用户数成千上万甚至过亿,存储用户属性的二维矩阵会占用巨大的存储空间,而商品数却没有那么多,所以这个模型的好处显而易见。但是它有个缺点,就是迭代时间太长了,这是可以预见的,以时间换空间嘛。

同样的,需要计算其各个参数的偏导数,求出更新式,显然,bu和bi的更新式和刚才求出的一样

其中的向量z等于下面这坨

现在要求qik、x和y的更新式,在这里如果把z看成Pu则根据前面求得的更新式可得qik的更新式:

求Xj的导数需要有点耐心,SSE的第二项对Xj的导数很容易得到,现在来求第一项对Xj的导数:

上式求和符中的i,j都是属于R(u)的,可以看到Zm只有当m==k时才与Xjk有关,所以

因为

所以

所以得到SSE对Xjk的导数为

同理可得SSE对Yjk的导数为

所以得到Xjk和Yjk的更新方程

这就叫ASVD。。。。。。

SVDPP

最后这个模型也是Koren文章中提到的,SVDPlusPlus(SVD++),它的预测式子为

这里的N(u)表示用户u行为记录(包括浏览的和评过分的商品集合)。看了ASVD的更新式推导过程再来看这个应该很简单,Puk和qik的更新式子不变,Yjk的更新式子和ASVD的更新式一样:

这些就是Koren在NetFlix大赛中用到的SVD算法,最后,还有一个需要提的:

对偶算法

将前面预测公式中的u和i调换位置得到其对偶算法,对于RSVD而言,u和i的位置是等价的、对称的,所以其对偶算法和其本身没有区别,但对于ASVD和SVD++则不同,有时候对偶算法得到的结果更加精确,并且,如果将对偶算法和原始算法的预测结果融合在一起的话,效果的提升会让你吃惊!

对偶的ASVD预测公式:

这里R(i)表示评论过商品i的用户集合,N(i)表示浏览过商品i但没有评论的用户集合。由于用户数量庞大,所以对偶的ASVD会占用很大空间,这里需要做取舍了。

对偶的SVD++预测公式:

这里N(i)表示对商品i有过行为(浏览或评分)的用户集合。

实现这一对偶的操作其实很简单,只要读取数据的时候把用户id和商品id对调位置即可,也就是将R矩阵转置后再训练。

暂时实现的算法就这些,代码都在github上了,有兴趣的可以看看,地址:https://github.com/jingchenUSTC/SVDRecommenderSystem

 

 

 

奇异值分解(Singular Value Decomposition)是线性代数中一种重要的矩阵分解,是矩阵分析中正规矩阵酉对角化的推广。在信号处理、统计学等领域有重要应用。奇异值分解则是谱分析理论在任意矩阵上的推广。 [1]
中文名 奇异值分解 外文名 Singular Value Decomposition 特 征 对矩阵的扰动不敏感 应 用矩阵分解 应用学科 线性代数 所属领域 信号处理、统计学
目录
1 基本介绍
2 理论描述
3 几何意义
4 应用
▪ 求伪逆
▪ 平行奇异值
▪ 矩阵近似值
5 编程
基本介绍编辑
奇异值分解在某些方面与对称矩阵或Hermite矩阵基于特征向量的对角化类似。然而这两种矩阵分解尽管有其相关性,但还是有明显的不同。谱分析的基础是对称阵特征向量的分解,而奇异值分解则是谱分析理论在任意矩阵上的推广。 [1]
理论描述编辑
假设M是一个m×n阶矩阵,其中的元素全部属于域 K,也就是实数域或复数域。如此则存在一个分解使得

其中U是m×m阶酉矩阵;Σ是半正定m×n阶对角矩阵;而V*,即V的共轭转置,是n×n阶酉矩阵。这样的分解就称作M的奇异值分解。Σ对角线上的元素Σi,其中i即为M的奇异值。
常见的做法是为了奇异值由大而小排列。如此Σ便能由M唯一确定了。(虽然U和V仍然不能确定)
直观的解释
在矩阵M的奇异值分解中
·U的列(columns)组成一套对M的正交”输入”或”分析”的基向量。这些向量是MM*的特征向量。
·V的列(columns)组成一套对M的正交”输出”的基向量。这些向量是M*M的特征向量。
·Σ对角线上的元素是奇异值,可视为是在输入与输出间进行的标量的”膨胀控制”。这些是M*M及MM*的奇异值,并与U和V的列向量相对应。
奇异值和奇异向量, 以及他们与奇异值分解的关系
一个非负实数σ是M的一个奇异值仅当存在Km 的单位向量u和Kn的单位向量v如下 :
其中向量u 和v分别为σ的左奇异向量和右奇异向量。
对于任意的奇异值分解,矩阵Σ的对角线上的元素等于M的奇异值。U和V的列分别是奇异值中的左、右奇异向量。因此,上述定理表明: [2]
(1)一个m × n的矩阵至多有 p = min(m,n)个不同的奇异值;
(2)总是可以找到在Km 的一个正交基U,组成M的左奇异向量;
(3)总是可以找到和Kn的一个正交基V,组成M的右奇异向量。
如果一个奇异值中可以找到两个左(或右)奇异向量是线性相关的,则称为退化。
非退化的奇异值具有唯一的左、右奇异向量,取决于所乘的单位相位因子eiφ(根据实际信号)。因此,如果M的所有奇异值都是非退化且非零,则它的奇异值分解是唯一的,因为U中的一列要乘以一个单位相位因子且同时V中相应的列也要乘以同一个相位因子。
根据定义,退化的奇异值具有不唯一的奇异向量。因为,如果u1和u2为奇异值σ的两个左奇异向量,则两个向量的任意规范线性组合也是奇异值σ一个左奇异向量,类似的,右奇异向量也具有相同的性质。因此,如果M 具有退化的奇异值,则它的奇异值分解是不唯一的。
几何意义编辑
因为U 和V 向量都是单位化的向量, 我们知道U的列向量u1,…,um组成了K空间的一组标准正交基。同样,V的列向量v1,…,vn也组成了K空间的一组标准正交基(根据向量空间的标准点积法则)。 [2]
线性变换T:即K → K,把向量Nx变换为Mx。考虑到这些标准正交基,这个变换描述起来就很简单了:T(vi) = σi ui, for i = 1,…,min(m,n),其中σi 是对角阵Σ中的第i个元素;当i > min(m,n)时,T(vi) = 0。
这样,SVD理论的几何意义就可以做如下的归纳:对于每一个线性映射T: K → K,T把K的第i个基向量映射为K的第i个基向量的非负倍数,然后将余下的基向量映射为零向量。对照这些基向量,映射T就可以表示为一个非负对角阵。
应用编辑
求伪逆
奇异值分解可以被用来计算矩阵的伪逆。若矩阵 M 的奇异值分解为 ,那么 M 的伪逆为 。 [2]
其中 是 的伪逆,并将其主对角线上每个非零元素都求倒数之后再转置得到的。求伪逆通常可以用来求解线性最小平方、最小二乘法问题。
平行奇异值
把频率选择性衰落信道进行分解。
矩阵近似值
奇异值分解在统计中的主要应用为主成分分析(PCA),一种数据分析方法,用来找出大量数据中所隐含的“模式”,它可以用在模式识别,数据压缩等方面。PCA算法的作用是把数据集映射到低维空间中去。 数据集的特征值(在SVD中用奇异值表征)按照重要性排列,降维的过程就是舍弃不重要的特征向量的过程,而剩下的特征向量组成的空间即为降维后的空间。

 

 

SVD++

基于潜在(隐藏)因子的推荐,常采用SVD或改进的SVD++

奇异值分解(SVD)

考虑CF中最为常见的用户给电影评分的场景,我们需要一个数学模型来模拟用户给电影打分的场景,比如对评分进行预测。

将评分矩阵U看作是两个矩阵的乘积:

其中,uxy 可以看作是user x对电影的隐藏特质y的热衷程度,而iyz可以看作是特质 y 在电影 z中的体现程度。那么上述模型的评分预测公式为:

q 和 p 分别对应了电影和用户在各个隐藏特质上的特征向量。

以上的模型中,用户和电影都体现得无差别,例如某些用户非常挑剔,总是给予很低的评分;或是某部电影拍得奇烂,恶评如潮。为了模拟以上的情况,需要引入 baseline predictor.

其中 μ 为所有评分基准,bi 为电影 i 的评分均值相对μ的偏移,bu 类似。注意,这些均为参数,需要通过训练得到具体数值,不过可以用相应的均值作为初始化时的估计。

模型参数bi,bu,qi,pu通过最优化下面这个目标函数获得:

可以用梯度下降方法或迭代的最小二乘算法求解。在迭代最小二乘算法中,首先固定pu优化qi,然后固定qi优化pu,交替更新。梯度下降方法中参数的更新式子如下(为了简便,把目标函数中的μ+bi+bu+q⊤ipu整体替换为r^ui):

其中α是更新步长。

SVD++

某个用户对某个电影进行了评分,那么说明他看过这部电影,那么这样的行为事实上蕴含了一定的信息,因此我们可以这样来理解问题:评分的行为从侧面反映了用户的喜好,可以将这样的反映通过隐式参数的形式体现在模型中,从而得到一个更为精细的模型,便是 SVD++.

其中 I(u) 为该用户所评价过的所有电影的集合,yj为隐藏的“评价了电影 j”反映出的个人喜好偏置。收缩因子取集合大小的根号是一个经验公式,并没有理论依据。

模型参数bi,bu,qi,pu,yj通过最优化下面这个目标函数获得:

与SVD方法类似,可以通过梯度下降算法进行求解。

NMF

非负矩阵分解

通常的矩阵分解会把一个大的矩阵分解为多个小的矩阵,但是这些矩阵的元素有正有负。而在现实世界中,比如图

像,文本等形成的矩阵中负数的存在是没有意义的,所以如果能把一个矩阵分解成全是非负元素是很有意义的。在

NMF中要求原始的矩阵的所有元素的均是非负的,那么矩阵可以分解为两个更小的非负矩阵的乘积,这个矩阵

有且仅有一个这样的分解,即满足存在性唯一性

 

Contents

 

   1. NMF问题描述

   2. NMF实现原理

   3. NMF的应用

   4. NMF的R实现

   5. NMF的Julia实现

   6. 结束语

 

 

1. NMF问题描述

 

传统的NMF问题可以描述如下

 

给定矩阵,寻找非负矩阵和非负矩阵,使得

 

分解前后可理解为:原始矩阵的列向量是对左矩阵中所有列向量的加权和,而权重系数就是右矩阵对应列

向量的元素,故称为基矩阵,为系数矩阵。一般情况下的选择要比小,即满足

这时用系数矩阵代替原始矩阵,就可以实现对原始矩阵进行降维,得到数据特征的降维矩阵,从而减少存储空

间,减少计算机资源。

 

 

2. NMF实现原理

 

NMF求解问题实际上是一个最优化问题,利用乘性迭代的方法求解,非负矩阵分解是一个NP问题。NMF

问题的目标函数有很多种,应用最广泛的就是欧几里得距离KL散度。

 

   在NMF的分解问题中,假设噪声矩阵为,那么有

 

 

现在要找出合适的使得最小。假设噪声服从不同的概率分布,通过最大似然函数会得到不同类型

的目标函数。接下来会分别以噪声服从高斯分布泊松分布来说明。

 

 

   (1)噪声服从高斯分布

 

假设噪声服从高斯分布,那么得到最大似然函数为

 

 

取对数后,得到对数似然函数为

 

 

假设各数据点噪声的方差一样,那么接下来要使得对数似然函数取值最大,只需要下面目标函数值最小。

 

 

该损失函数为2范数损失函数,它是基于欧几里得距离的度量。又因为

 

 

那么得到

 

 

同理有

 

 

接下来就可以使用梯度下降法进行迭代了。如下

 

 

如果选取

 

 

那么最终得到迭代式为

 

 

可看出这是乘性迭代规则,每一步都保证了结果为正数,一直迭代下去就会收敛,当然收敛性的证明省略。

 

 

   (2)噪声服从泊松分布

 

       若噪声为泊松噪声,那么得到损失函数为

 

 

同样经过推到得到

 

 

3. NMF的应用

 

NMF能用于发现数据库中图像的特征,便于快速识别应用,比如实现录入恐怖分子的照片,然后在安检口对可疑

人员进行盘查。在文档方面,NMF能够发现文档的语义相关度,用于信息的自动索引和提取。在生物学中,在

DNA阵列分析中识别基因等。在语音识别系统中NMF也能发挥重要作用。

 

 

4. NMF的R实现

 

   先说NMF的R使用,R语言中已经有NMF包可用于NMF。非负矩阵分解已经有了很多算法,例如Multiplicative

Update,Projected Gradient Method,Gradient Descent。交替最小二乘法比较符合上面的解释,具

有更好的统计意义,而且收敛性质很好,计算速度快。已有的交替最小二乘法都使用冷冰冰的二次规划算法求解

非负回归,在bignmf中使用的最小角回归的思路解非负回归,更有统计的感觉,而且速度更快。安装包如下

 

bignmf:https://github.com/panlanfeng/bignmf

 

bignmf的部分源码如下

 

 

从上面代码可以看出,当矩阵V中元素是非double时,会强制转化为double类型。而WH矩阵是服从高斯分布

的随机矩阵,随后传入参数调用whupdate,如果迭代次数太少会打印警告Iteration doesn’t converge!

 

所以用bignmf分解非负矩阵调用如下函数就行了

 

 

代码:

 

以上是NMF分解的R语言实现,bignmf借助伟大的Rcpp和RcppEigen包,算法内层用C++ 代码实现接下来会接

着讲NMF在Julia中的实现。

 

 

5. NMF的Julia实现

 

在Julia语言中,有一个NMF包专门用来进行NMF分解的,现在就来详细了解它。首先来安装NMF,在Julia交互

式环境下使用如下命令

 

 

在Julia中,关于NMF的计算有很多方法,详见:https://github.com/JuliaStats/NMF.jl

 

NMF分解的一个例子如下

 

 

 

6. 结束语

 

关于NMF的解的收敛性和唯一性以后再证明,除了NMF分解以外,还有一种非负分解也用的比较多,而且很有

效,叫做负张量分解法。可以参考如下论文

 

 

Slope One

推荐算法之 slope one 算法

1.示例引入

多个吃货在某美团的某家饭馆点餐,如下两道菜:

可乐鸡翅:
这里写图片描述
红烧肉:
这里写图片描述

顾客吃过后,会有相关的星级评分。假设评分如下:
评分 可乐鸡翅 红烧肉
小明 4 5
小红 4 3
小伟 2 3
小芳 3 ?
问题:请猜测一下小芳可能会给“红烧肉”打多少分?
思路:把两道菜的平均差值求出来,可乐鸡翅减去红烧肉的平均偏差:[(4-5)+(4-3)+(2-3)]/3=-0.333。一个新客户比如小芳,只吃了可乐鸡翅评分为3分,那么可以猜测她对红烧肉的评分为:3-(-0.333)=3.333

这就是slope one 算法的基本思路,非常非常的简单。

2.slope one 算法思想

Slope One 算法是由 Daniel Lemire 教授在 2005 年提出的一个Item-Based 的协同过滤推荐算法。和其它类似算法相比, 它的最大优点在于算法很简单, 易于实现, 执行效率高, 同时推荐的准确性相对较高。
Slope One算法是基于不同物品之间的评分差的线性算法,预测用户对物品评分的个性化算法。主要两步:
Step1:计算物品之间的评分差的均值,记为物品间的评分偏差(两物品同时被评分);
这里写图片描述
Step2:根据物品间的评分偏差和用户的历史评分,预测用户对未评分的物品的评分。
这里写图片描述
Step3:将预测评分排序,取topN对应的物品推荐给用户。

举例:
假设有100个人对物品A和物品B打分了,R(AB)表示这100个人对A和B打分的平均偏差;有1000个人对物品B和物品C打分了, R(CB)表示这1000个人对C和B打分的平均偏差;
这里写图片描述

3.python实现

3.1数据

def loadData():
    items={'A':{1:5,2:3},
           'B':{1:3,2:4,3:2},
           'C':{1:2,3:5}}
    users={1:{'A':5,'B':3,'C':2},
           2:{'A':3,'B':4},
           3:{'B':2,'C':5}}
    return items,users

3.2物品间评分偏差

#***计算物品之间的评分差
#items:从物品角度,考虑评分
#users:从用户角度,考虑评分
def buildAverageDiffs(items,users,averages):
    #遍历每条物品-用户评分数据
    for itemId in items:
        for otherItemId in items:
            average=0.0 #物品间的评分偏差均值
            userRatingPairCount=0 #两件物品均评过分的用户数
            if itemId!=otherItemId: #若无不同的物品项
                for userId in users: #遍历用户-物品评分数
                    userRatings=users[userId] #每条数据为用户对物品的评分
                    #当前物品项在用户的评分数据中,且用户也对其他物品由评分
                    if itemId in userRatings and otherItemId in userRatings:
                        #两件物品均评过分的用户数加1
                        userRatingPairCount+=1
                        #评分偏差为每项当前物品评分-其他物品评分求和
                        average+=(userRatings[otherItemId]-userRatings[itemId])
                averages[(itemId,otherItemId)]=average/userRatingPairCount

3.3预估评分

#***预测评分
#users:用户对物品的评分数据
#items:物品由哪些用户评分的数据
#averages:计算的评分偏差
#targetUserId:被推荐的用户
#targetItemId:被推荐的物品
def suggestedRating(users,items,averages,targetUserId,targetItemId):
    runningRatingCount=0 #预测评分的分母
    weightedRatingTotal=0.0 #分子
    for i in users[targetUserId]:
        #物品i和物品targetItemId共同评分的用户数
        ratingCount=userWhoRatedBoth(users,i,targetItemId)
        #分子
        weightedRatingTotal+=(users[targetUserId][i]-averages[(targetItemId,i)])\
        *ratingCount
        #分母
        runningRatingCount+=ratingCount
    #返回预测评分
    return weightedRatingTotal/runningRatingCount

统计两物品共同评分的用户数

# 物品itemId1与itemId2共同有多少用户评分
def userWhoRatedBoth(users,itemId1,itemId2):
    count=0
    #用户-物品评分数据
    for userId in users:
        #用户对物品itemId1与itemId2都评过分则计数加1
        if itemId1 in users[userId] and itemId2 in users[userId]:
            count+=1
    return count

3.4测试结果:

if __name__=='__main__':
    items,users=loadData()
    averages={}
    #计算物品之间的评分差
    buildAverageDiffs(items,users,averages)
    #预测评分:用户2对物品C的评分
    predictRating=suggestedRating(users,items,averages,2,'C')
    print 'Guess the user will rate the score :',predictRating

结果:用户2对物品C的预测分值为
Guess the user will rate the score : 3.33333333333

4.slopeOne使用场景

该算法适用于物品更新不频繁,数量相对较稳定并且物品数目明显小于用户数的场景。依赖用户的用户行为日志和物品偏好的相关内容。
优点:
1.算法简单,易于实现,执行效率高;
2.可以发现用户潜在的兴趣爱好;
缺点:
依赖用户行为,存在冷启动问题和稀疏性问题。

 

 

k-NN
Centered k-NN
k-NN Baseline

 

Co-Clustering
Baseline

 

 

基于协同过滤,NMF和Baseline的推荐算法

转载 2016年04月05日 04:59:52

 老早就想整理一篇推荐算法的入门博文,今天抽空写一下。本文以电影推荐系统为例,简单地介绍基于协同过滤,PMF概率矩阵分解,NMF非负矩阵分解和Baseline的推荐系统算法。NMF的实现具体可以参考Reference中的「基于矩阵分解的推荐算法,简单入门」一文,对我启发很大。

杂谈

老早就想整理一篇推荐算法的入门博文,今天抽空写一下。本文以电影推荐系统为例,简单地介绍基于协同过滤,PMF概率矩阵分解,NMF非负矩阵分解和Baseline的推荐系统算法。NMF的实现具体可以参考Reference中的「基于矩阵分解的推荐算法,简单入门」一文,对我启发很大。

基于协同过滤的推荐算法

什么是协同过滤?协同过滤是利用集体智慧的一个典型方法。要理解什么是协同过滤 (Collaborative Filtering, 简称 CF),首先想一个简单的问题,如果你现在想看个电影,但你不知道具体看哪部,你会怎么做?大部分的人会问问周围的朋友,看看最近有什么好看的电影推荐,而我们一般更倾向于从口味比较类似的朋友那里得到推荐;或者,搜索与你喜欢的电影同类型的电影推荐。

User-based的推荐算法

如上图,我们收集到用户-电影评价矩阵,假设用户A对于物品D的评价为null,这时我们对比用户A、用户B、用户C的特征向量(以物品评价为特征),可以发现用户A与用户C的相似度较大,这时我们可以认为,对于用户C喜欢的物品D,用户A也应该喜欢它,这是就把物品D推荐给用户A。

Item-based的推荐算法

同理,我们对比物品A、物品B、物品C的特征向量(以用户对该物品的喜欢程度为特征),发现物品A与物品C很像,就把用品C推荐给喜欢物品A的用户C。

SVD和PMF概率矩阵分解

在实际业务场景中,user-Item矩阵有可能非常稀疏,存储率有可能连1%都达不到。怎么办呢?通常使用矩阵分解算法来提取出更有用的信息。SVD在矩阵分解方面可以参考基于SVD实现PCA的图像识别 一文,把它分解成用户矩阵(左奇异特征向量矩阵)和物品矩阵(右奇异特征向量矩阵)分别代表各自的特性。

但是SVD算法的时间复杂度很大,不适合用来解决这种比数据量较大的问题,这时就有PMF概率矩阵分解。把用户-电影 评分看成一个矩阵,rui表示u对电影i的评分,于是电影评分矩阵可以这样来估计:

其中P和Q就相当于SVD中的前k 个特征向量构成的矩阵,分别描述user-based和item-based。在PMF中使用SGD(随机梯度下降)进行优化时,使用如下的迭代公式:

我们把证明放到下一章节 NMF非负矩阵分解,NMF其实就是在PMF的基础上加入一点约束,具体约束公式如下:

基于NMF非负矩阵分解的推荐算法

我们知道,要做推荐系统,最基本的一个数据就是,用户-物品的评分矩阵,如下图所示:

矩阵中,描述了5个用户(U1,U2,U3,U4 ,U5)对4个物品(D1,D2,D3,D4)的评分(1-5分),“-” 表示没有评分,现在目的是把没有评分的 给预测出来,然后按预测的分数高低,给用户进行推荐。

如何预测缺失的评分呢?对于缺失的评分,可以转化为基于机器学习的回归问题,也就是连续值的预测,对于矩阵分解有如下式子,R是类似图的评分矩阵,假设N*M维(N表示行数,M表示列数),可以分解为P跟Q矩阵,其中P矩阵维度N*K,P矩阵维度K*M。

对于P,Q矩阵的解释,直观上,P矩阵是N个用户对K个主题的关系,Q矩阵是K个主题跟M个物品的关系,至于K个主题具体是什么,在算法里面K是一个参数,需要调节的,通常10~100之间。

对于上式的左边项,表示的是R^ 第i 行,第j 列的元素值,对于如何衡量矩阵分解的好坏,我们给出如下风险函数:

有了风险函数,我们就可以采用梯度下降法不断地减小损失值,直到不能再减小为止,最后的目标,就是每一个元素(非缺失值)的e(i,j)的总和 最小。我们可以得到如下梯度以及p、q的更新方式(其中α是学习步长,详见李航《统计机器学习》):

在训练p、q参数过程中,为了防止过拟合,我们给出一个正则项,风险函数修改如下:

相应的p、q参数学习更新方式如下:

至此,我们就可以学习出p、q矩阵,将p x q就可以得到新的估计矩阵,由于加入了非负处理(缺失值部分的处理),我们可以发现原先缺失的地方有了一个估计值,这个估计值就作为了推荐的分值(其实就是拿非缺失值部分作参数学习训练,学习出来的结果当然不会有负数)。NMF实现代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
packagenmf;
publicclassNmf {
    publicdouble[][] RM, PM, QM;
    publicintKc, Uc, Oc;
    publicintsteps;
    publicdoubleAlpha, Beta;
    publicvoidrun() {
        for(ints = 0; s < steps; s++) {
            // 梯度下降更新
            for(inti = 0; i < Uc; i++) {
                for(intj = 0; j < Oc; j++) {
                    if(RM[i][j] > 0) {
                        // 计算eij
                        doublee = 0, pq = 0;
                        for(intk = 0; k < Kc; k++) {
                            pq += PM[i][k] * QM[k][j];
                        }
                        e = RM[i][j] - pq;
                        // 更新Pik和Qkj,同时保证非负
                        for(intk = 0; k < Kc; k++) {
                            PM[i][k] += Alpha
                                    * (2* e * QM[k][j] - Beta * PM[i][k]);
                            // PM[i][k] = PM[i][k] > 0 ? PM[i][k] : 0;
                            QM[k][j] += Alpha
                                    * (2* e * PM[i][k] - Beta * QM[k][j]);
                            // QM[k][j] = QM[k][j] > 0 ? QM[k][j] : 0;
                        }
                    }
                }
            }
            // 计算风险损失
            doubleloss = 0;
            for(inti = 0; i < Uc; i++) {
                for(intj = 0; j < Oc; j++) {
                    if(RM[i][j] > 0) {
                        // 计算eij^2
                        doublee2 = 0, pq = 0;
                        for(intk = 0; k < Kc; k++) {
                            pq += PM[i][k] * QM[k][j];
                        }
                        e2 = Math.pow(RM[i][j] - pq, 2);
                        for(intk = 0; k < Kc; k++) {
                            e2 += Beta
                                    2
                                    * (Math.pow(PM[i][k], 2) + Math.pow(
                                            QM[k][j], 2));
                        }
                        loss += e2;
                    }
                }
            }
            if(loss < 0.01) {
                System.out.println("OK");
                break;
            }
            // if (s % 100 == 0) {
            // System.out.println(loss);
            // }
        }
    }
    publicNmf(double[][] RM, double[][] PM, double[][] QM, intKc, intUc,
            intOc, intsteps, doubleAlpha, doubleBeta) {
        this.RM = RM;
        this.PM = PM;
        this.QM = QM;
        this.Kc = Kc;
        this.Uc = Uc;
        this.Oc = Oc;
        this.steps = steps;
        this.Alpha = Alpha;
        this.Beta = Beta;
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
packagenmf;
importjava.util.Scanner;
publicclassKeyven {
    publicstaticvoidmain(String[] args) {
        intUc = 5, Oc = 4, Kc = 2;
        double[][] RM = newdouble[Uc][Oc];
        double[][] PM = newdouble[Uc][Kc];
        double[][] QM = newdouble[Kc][Oc];
        /*
         * 5 3 0 1 4 0 0 1 1 1 0 5 1 0 0 4 0 1 5 4
         */
        Scanner input = newScanner(System.in);
        for(inti = 0; i < Uc; i++) {
            for(intj = 0; j < Oc; j++) {
                RM[i][j] = input.nextDouble();
            }
        }
        
        for(inti = 0; i < Uc; i++) {
            for(intj = 0; j < Oc; j++) {
                System.out.printf("%.2f\t", RM[i][j]);
            }
            System.out.println();
        }
        System.out.println();
        for(inti = 0; i < Uc; i++) {
            for(intj = 0; j < Kc; j++) {
                PM[i][j] = Math.random() % 9;
            }
        }
        for(inti = 0; i < Kc; i++) {
            for(intj = 0; j < Oc; j++) {
                QM[i][j] = Math.random() % 9;
            }
        }
        // 最多迭代5000次,学习步长控制为0.002,正则项参数设置为0.02
        Nmf nmf = newNmf(RM, PM, QM, Kc, Uc, Oc, 50000.0020.02);
        nmf.run();
        for(inti = 0; i < Uc; i++) {
            for(intj = 0; j < Oc; j++) {
                doubletemp = 0;
                for(intk = 0; k < Kc; k++) {
                    temp += PM[i][k] * QM[k][j];
                }
                System.out.printf("%.2f\t", temp);
            }
            System.out.println();
        }
        input.close();
    }
}

实验结果:

基于Baseline的推荐算法

要评估一个策略的好坏,就需要建立一个对比基线,以便后续观察算法效果的提升。此处我们可以简单地对推荐算法进行建模作为基线。假设我们的训练数据为:  <user, item, rating>三元组, 其中user为用户id, item为物品id(item可以是MovieLens上的电影,Amazon上的书, 或是百度关键词工具上的关键词), rating为user对item的投票分数, 其中用户u对物品i的真实投票分数我们记为rui,基线(baseline)模型预估分数为bui,则可建模如下:

其中mu(希腊字母mu)为所有已知投票数据中投票的均值,bu为用户的打分相对于平均值的偏差(如果某用户比较苛刻,打分都相对偏低, 则bu会为负值;相反,如果某用户经常对很多片都打正分, 则bu为正值); bi为该item被打分时,相对于平均值得偏差,可反映电影受欢迎程度。 bui则为基线模型对用户u给物品i打分的预估值。该模型虽然简单, 但其中其实已经包含了用户个性化和item的个性化信息, 而且特别简单(很多时候, 简单就是一个非常大的特点, 特别是面对大规模数据时)。

基线模型中, mu可以直接统计得到,我们的优化函数可以写为(其实就是最小二乘法):

也可以直接写成如下式子,因为它本身就是经验似然:

上述式子中u∈R(i) 表示评价过电影 i 的所有用户,|R(i)| 为其集合的个数;同理,i∈R(u) 表示用户 u 评价过的所有电影,|R(u)| 为其集合的个数。实现代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
packagebaseline;
publicclassBaseline {
    publicdouble[] bi, bu;
    publicdouble[][] RM;
    publicintUc, Ic;
    publicdoublelamada2, lamada3;
    publicBaseline(double[][] RM, intUc, intIc, doublelamada2,
            doublelamada3) {
        this.RM = RM;
        this.lamada2 = lamada2;
        this.lamada3 = lamada3;
        this.Uc = Uc;
        this.Ic = Ic;
        this.bu = newdouble[Uc];
        this.bi = newdouble[Ic];
    }
    publicvoidrun() {
        // 计算μ
        doubleavg = 0;
        for(inti = 0; i < Uc; i++) {
            for(intj = 0; j < Ic; j++) {
                avg += RM[i][j];
            }
        }
        avg = avg / Uc / Ic;
        // 更新bi
        for(inti = 0; i < Ic; i++) {
            doublebis = 0;
            intIcnt = 0// 点评过电影i的所有User个数
            for(inttu = 0; tu < Uc; tu++) {
                if(RM[tu][i] != 0) {
                    bis += RM[tu][i] - avg;
                    Icnt++;
                }
            }
            bi[i] = bis / ((double)Icnt + lamada2);
        }
        // 更新bu
        for(intu = 0; u < Uc; u++) {
            doublebus = 0;
            intUcnt = 0// 用户u点评过得电影Item个数
            for(intti = 0; ti < Ic; ti++) {
                if(RM[u][ti] != 0) {
                    bus += RM[u][ti] - avg - bi[ti];
                    Ucnt++;
                }
            }
            bu[u] = bus / ((double)Ucnt + lamada3);
        }
        for(intu = 0; u < Uc; u++) {
            for(inti = 0; i < Ic; i++) {
                if(RM[u][i] == 0) {
                    RM[u][i] = avg + bi[i] + bu[u];
                }
            }
        }
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
packagebaseline;
importjava.util.Scanner;
publicclassKeyven {
    publicstaticvoidmain(String[] args) {
        intUc = 5, Ic = 4;
        double[][] RM = newdouble[Uc][Ic];
        /*
         * 5 3 0 1 4 0 0 1 1 1 0 5 1 0 0 4 0 1 5 4
         */
        Scanner input = newScanner(System.in);
        for(inti = 0; i < Uc; i++) {
            for(intj = 0; j < Ic; j++) {
                RM[i][j] = input.nextDouble();
            }
        }
        Baseline bl = newBaseline(RM, Uc, Ic, 00);
        bl.run();
        for(inti = 0; i < Uc; i++) {
            for(intj = 0; j < Ic; j++) {
                System.out.printf("%.2f\t", RM[i][j]);
            }
            System.out.println();
        }
        input.close();
    }
}

Baseline基线模型与NMF矩阵分解模型试验效果对比如下:

后记

什么是梯度下降?考虑下图一种简单的情形,风险函数为loss = kx,则总体损失就是积分∫kx dx,取梯度的反方向进行逐步更新至总体损失减小… …在我看来,其实,数据挖掘 = ①线性代数+②应用概率统计+③高数(积分、梯度等数理意义)+④李航《统计机器学习》神书+⑤算法与数据结构… …只要好好努力打好基础,人就能不断向前走下去^_^

Reference

从item-base到svd再到rbm,多种Collaborative Filtering(协同过滤算法)从原理到实现 

百度电影推荐系统比赛 小结 ——记我的初步推荐算法实践

白话NMF(Non-negative Matrix Factorization)

基于矩阵分解的推荐算法,简单入门(证明算法,非常有用!)

SVD因式分解实现协同过滤-及源码实现

探索推荐引擎内部的秘密,第 2 部分: 深入推荐引擎相关算法 – 协同过滤 

推荐系统中近邻算法与矩阵分解算法效果的比较——基于movielens数据集

 

 

 

 

 

Random

Leave a Reply

Your email address will not be published. Required fields are marked *