机器学习实战及Python实现——奇异值分解（SVD）实现简单推荐系统

本篇讲数据降维的另一种更普遍的算法——奇异值分解，主要内容包括数学原理，计算步骤，优缺点，应用场景、Python推荐示例等内容。

1、数学原理

奇异值分解（Singular Value Decomposition，SVD）是一种重要的矩阵分解。与之相对的是特征值分解（主成分分析主要使用方法），但特征值分解是针对的是方阵，但在实际应用场景中，我们经常遇到的矩阵都不是方阵，比如N个学生，每个学生的M科成绩，其中N≠M，这就组成N*M的非方阵矩阵。

对于一般普通的矩阵（包括方阵矩阵），如何来描述其重要特征？奇异值分解就是来做这些事情的。其中的矩阵分解公式为：

假设A是一个M* N的矩阵，那么通过矩阵分解将会得到U，Σ，V’（V的转置）三个矩阵，其中U是一个M * M的方阵，被称为左奇异向量，方阵里面的向量是正交的；Σ是一个M* N的对角矩阵，除了对角线的元素其他都是0，对角线上的值称为奇异值；V’(V的转置)是一个N * N的矩阵，被称为右奇异向量，方阵里面的向量也都是正交的。用图形展示如下图：

2、计算步骤

该矩阵是如何分解的？奇异值和特征值是如何计算的？

（1）将矩阵A的转置 * A，将会得到一个方阵，将方阵进行特征值分解：

其中得到的v，就是右奇异向量。

（2）通过方阵还可以求解σ和u：

σ是上文提到的奇异值，u是上文提到的左奇异向量。其中奇异值σ跟特征值很类似，在矩阵Σ中也是从大到小排列，而且σ的减少特别的快，在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上了。也就是说，我们也可以用前r大的奇异值来近似描述矩阵，其中r<<n，这里定义一下部分奇异值分解：

（3）选择适当的r，其中r是一个远小于m、n的数，这样可将原矩阵分解为：

其中，右边的三个矩阵相乘的结果将会是一个接近于A的矩阵，而r越接近n，其相乘结果越接近A。根据储存原理，储存量与矩阵面积正相关，因此面积越小占用的储存空间越小。而三个矩阵的面积之和要远小于原矩阵。因此如果要储存A的信息，只需要储存U、Σ、V就可以，因此信息得到压缩。

3、优缺点

奇异值分解的优点是：可以简化数据，压缩维度，去除数据噪音，提升算法的结果，加快模型计算性能，可以针对任一普通矩阵进行分解（包括样本数小于特征数），不受限于方阵。

奇异值分解的缺点是：转换后的数据比较难理解，如何与具体业务知识对应起来是难点。

4、应用场景

奇异值分解应用场景一：隐性语义索引（Latent Semantic Indexing，LSI）

矩阵是有文档（M行）和词语（N列）组成，通过奇异值分解，可以分析出那些文档或词语属于同一主题或概念，可应用于更高效的文档检索

奇异值分解应用场景二：推荐系统

通过奇异值分解，可以计算项与人之间的相似度，而进行协同过滤，向用户推荐相关产品。

5、数据背景

（1）样本数据

本次样本数据是11*11，其中行表示用户，列表示食品，中间数字表示该用户对食品的打分。如果数字为0，表示该用户没有吃过该食品。本次模型的目的就是向用户推荐未吃过的食品。

（2）推荐思路

首先，寻找用户未评价的食品，即用户-矩阵中的0值；

再次，对用户未打分的食品，通过相似度计算预计其可能会打多少分数；

最后，对这些打分的食品根据评分从高到低进行排序，返回前N个食品，这就是推荐结果。

（3）相似度计算

如何来衡量两个物品之间的相似情况，一般有以下三种方法

第一种是：欧氏距离

示例：

为将距离映射到【0,1】中，相似度=1/(1+欧氏距离)

第二种是：皮尔森相关系数

示例：

皮尔森系数在【-1,1】之间，为映射到【0,1】之间，相似度=0.5+0.5*corroef

第三种是：余弦夹角

余弦夹角在【-1,1】之间，为映射到【0,1】之间，相似度=0.5+0.5*cos

6、具体Python实现

1. from numpy import *  
2. from numpy import linalg as la  
3.   
4. def eulidSim(inA,inB):  
5.     return 1.0/(1.0+la.norm(inA,inB))  
6.   
7. def pearsSim(inA,inB):  
8.     if len(inA<3):return 1.0  
9.     return 0.5+0.5*corrcoef(inA,inB,rowvar=0)[0][1]  
10.   
11. def cosSim(inA,inB):  
12.     num=float(inA.T*inB)  
13.     denom=la.norm(inA)*la.norm(inB)  
14.     return 0.5+0.5*(num/denom)  

1. from numpy import *  
2. from numpy import linalg as la  
3.   
4. def loadExData():  
5.   return[[1, 1, 1, 0, 0],  
6.     [2, 2, 2, 0, 0],  
7.     [1, 1, 1, 0, 0],  
8.     [5, 5, 5, 0, 0],  
9.     [1, 1, 0, 2, 2],  
10.     [0, 0, 0, 3, 3],  
11.     [0, 0, 0, 1, 1]]  
12.       
13. def loadExData2():  
14.     return[[0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 5],  
15.            [0, 0, 0, 3, 0, 4, 0, 0, 0, 0, 3],  
16.            [0, 0, 0, 0, 4, 0, 0, 1, 0, 4, 0],  
17.            [3, 3, 4, 0, 0, 0, 0, 2, 2, 0, 0],  
18.            [5, 4, 5, 0, 0, 0, 0, 5, 5, 0, 0],  
19.            [0, 0, 0, 0, 5, 0, 1, 0, 0, 5, 0],  
20.            [4, 3, 4, 0, 0, 0, 0, 5, 5, 0, 1],  
21.            [0, 0, 0, 4, 0, 4, 0, 0, 0, 0, 4],  
22.            [0, 0, 0, 2, 0, 2, 5, 0, 0, 1, 2],  
23.            [0, 0, 0, 0, 5, 0, 0, 0, 0, 4, 0],  
24.            [1, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0]]  
25.       
26. def ecludSim(inA,inB):  
27.     return 1.0/(1.0 + la.norm(inA - inB))  
28.   
29. def pearsSim(inA,inB):  
30.     if len(inA) < 3 : return 1.0  
31.     return 0.5+0.5*corrcoef(inA, inB, rowvar = 0)[0][1]  
32.   
33. def cosSim(inA,inB):  
34.     num = float(inA.T*inB)  
35.     denom = la.norm(inA)*la.norm(inB)  
36.     return 0.5+0.5*(num/denom)  
37.   
38.   
39. #计算在给定相似度计算方法的条件下，用户对物品的估计评分值  
40. #standEst()函数中：参数dataMat表示数据矩阵，user表示用户编号，simMeas表示相似度计算方法，item表示物品编号  
41. def standEst(dataMat,user,simMeas,item):  
42.     n=shape(dataMat)[1] #shape用于求矩阵的行列  
43.     simTotal=0.0; ratSimTotal=0.0  
44.     for j in range(n):  
45.         userRating=dataMat[user,j]  
46.         if userRating==0:continue #若某个物品评分值为0，表示用户未对物品评分，则跳过，继续遍历下一个物品  
47.         #寻找两个用户都评分的物品  
48.         overLap=nonzero(logical_and(dataMat[:,item].A>0,dataMat[:,j].A>0))[0]  
49.   
50.         if len(overLap)==0:similarity=0  
51.         else: similarity=simMeas(dataMat[overLap,item],dataMat[overLap,j])  
52.   
53.         #print'the %d and%d similarity is: %f' %(item,j,similarity)  
54.         simTotal+=similarity  
55.         ratSimTotal+=similarity*userRating  
56.     if simTotal==0: return 0  
57.     else: return ratSimTotal/simTotal  
58.   
59. def recommend(dataMat,user,N=3,simMeas=cosSim,estMethod=standEst):  
60.     #寻找未评级的物品  
61.     unratedItems=nonzero(dataMat[user,:].A==0)[1]  
62.     if len(unratedItems)==0: return 'you rated everything'  
63.     itemScores=[]  
64.     for item in unratedItems:  
65.         estimatedScore=estMethod(dataMat,user,simMeas,item) #对每一个未评分物品，调用standEst()来产生该物品的预测得分  
66.         itemScores.append((item,estimatedScore)) #该物品的编号和估计得分值放入一个元素列表itemScores中  
67.     #对itemScores进行从大到小排序，返回前N个未评分物品  
68.     return sorted(itemScores,key=lambda jj:jj[1],reverse=True)[:N]  
69.   
70. def svdEst(dataMat, user, simMeas, item):  
71.     n = shape(dataMat)[1]  
72.     simTotal = 0.0; ratSimTotal = 0.0  
73.     U,Sigma,VT = la.svd(dataMat)  
74.     Sig4 = mat(eye(4)*Sigma[:4]) #arrange Sig4 into a diagonal matrix  
75.     xformedItems = dataMat.T * U[:,:4] * Sig4.I  #create transformed items  
76.     for j in range(n):  
77.         userRating = dataMat[user,j]  
78.         if userRating == 0 or j==item: continue  
79.         similarity = simMeas(xformedItems[item,:].T,\  
80.                              xformedItems[j,:].T)  
81.         print 'the %d and %d similarity is: %f' % (item, j, similarity)  
82.         simTotal += similarity  
83.         ratSimTotal += similarity * userRating  
84.     if simTotal == 0: return 0  
85.     else: return ratSimTotal/simTotal  

其中dataMat[:,item].A，表示找出item列，因为是matrix，用.A转成array，logical_and，其实就是找出最item列和j列都>0，只有都大于0才会是true，nonzero会给出其中不为0的index。

进行SVD分解：

[python]  view plain  copy

1. >>>from numpy import linalg as la  
2. >>> U,Sigma,VT=la.svd(mat(svdRec.loadExData2()))  
3. >>> Sigma  
4. array([ 1.38487021e+01, 1.15944583e+01, 1.10219767e+01,  
5.         5.31737732e+00, 4.55477815e+00, 2.69935136e+00,  
6.         1.53799905e+00, 6.46087828e-01, 4.45444850e-01,  
7.         9.86019201e-02, 9.96558169e-17])  

如何决定r？有个定量的方法是看多少个奇异值可以达到90%的能量，其实和PCA一样，由于奇异值其实是等于data×dataT特征值的平方根，所以总能量就是特征值的和

[python]  view plain  copy

1. >>> Sig2=Sigma**2  
2. >>> sum(Sig2)  
3. 541.99999999999932  

而取到前4个时，发现总能量大于90%，因此r=4

[python]  view plain  copy

1. >>> sum(Sig2[:3])  
2. 500.50028912757909  

SVD分解的关键在于，降低了user的维度，从n变到了4

[python]  view plain  copy

1. def svdEst(dataMat, user, simMeas, item):  
2.     n = shape(dataMat)[1]  
3.     simTotal = 0.0; ratSimTotal = 0.0  
4.     U,Sigma,VT = la.svd(dataMat)  
5.     Sig4 = mat(eye(4)*Sigma[:4]) #arrange Sig4 into a diagonal matrix  
6.     xformedItems = dataMat.T * U[:,:4] * Sig4.I  #create transformed items  
7.     for j in range(n):  
8.         userRating = dataMat[user,j]  
9.         if userRating == 0 or j==item: continue  
10.         similarity = simMeas(xformedItems[item,:].T,\  
11.                              xformedItems[j,:].T)  
12.         print 'the %d and %d similarity is: %f' % (item, j, similarity)  
13.         simTotal += similarity  
14.         ratSimTotal += similarity * userRating  
15.     if simTotal == 0: return 0  
16.     else: return ratSimTotal/simTotal  

其中关键一步，dataMat.T * U[:,:4] * Sig4.I

将m×n的dataMat用特征值缩放转换为n×4的item和user类的矩阵

[python]  view plain  copy

1. >>> myMat=mat(svdRec.loadExData2())  
2. >>> myMat  
3. matrix([[0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 5],  
4.         [0, 0, 0, 3, 0, 4, 0, 0, 0, 0, 3],  
5.         [0, 0, 0, 0, 4, 0, 0, 1, 0, 4, 0],  
6.         [3, 3, 4, 0, 0, 0, 0, 2, 2, 0, 0],  
7.         [5, 4, 5, 0, 0, 0, 0, 5, 5, 0, 0],  
8.         [0, 0, 0, 0, 5, 0, 1, 0, 0, 5, 0],  
9.         [4, 3, 4, 0, 0, 0, 0, 5, 5, 0, 1],  
10.         [0, 0, 0, 4, 0, 4, 0, 0, 0, 0, 4],  
11.         [0, 0, 0, 2, 0, 2, 5, 0, 0, 1, 2],  
12.         [0, 0, 0, 0, 5, 0, 0, 0, 0, 4, 0],  
13.         [1, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0]])  
14. >>> svdRec.recommend(myMat,1,estMethod=svdRec.svdEst)  
15. the 0 and 3 similarity is: 0.490950  
16. the 0 and 5 similarity is: 0.484274  
17. the 0 and 10 similarity is: 0.512755  
18. the 1 and 3 similarity is: 0.491294  
19. the 1 and 5 similarity is: 0.481516  
20. the 1 and 10 similarity is: 0.509709  
21. the 2 and 3 similarity is: 0.491573  
22. the 2 and 5 similarity is: 0.482346  
23. the 2 and 10 similarity is: 0.510584  
24. the 4 and 3 similarity is: 0.450495  
25. the 4 and 5 similarity is: 0.506795  
26. the 4 and 10 similarity is: 0.512896  
27. the 6 and 3 similarity is: 0.743699  
28. the 6 and 5 similarity is: 0.468366  
29. the 6 and 10 similarity is: 0.439465  
30. the 7 and 3 similarity is: 0.482175  
31. the 7 and 5 similarity is: 0.494716  
32. the 7 and 10 similarity is: 0.524970  
33. the 8 and 3 similarity is: 0.491307  
34. the 8 and 5 similarity is: 0.491228  
35. the 8 and 10 similarity is: 0.520290  
36. the 9 and 3 similarity is: 0.522379  
37. the 9 and 5 similarity is: 0.496130  
38. the 9 and 10 similarity is: 0.493617  
39. [(4, 3.3447149384692283), (7, 3.3294020724526967), (9, 3.328100876390069)]