机器学习笔记（Chapter 14 - SVD简化）

从数据中提取一些关键信息可以使用奇异值分解（Singular Value Decomposition，SVD），可以简化数据，去除噪声，将数据映射到低维空间。

SVD应用

• 奇异值分解的优点是简化数据，去除噪声，提高算法结果，用小得多的数据集表示原始数据集，实际上是取出了噪声。缺点是数据的转换可能难以理解。适用于数值型数据。
• 隐性语义索引（Latent Semantic Indexing，LSI）是SVD最早的应用之一。LSI中，一个矩阵是由文档和词语组成的，在该矩阵上应用SVD时，会构建出多个奇异值。这些奇异值代表了文档中的概念或主题，该特点可以用于更高效的文档搜索。在词语拼写错误或者出现同义词时，只基于词语存在与否的搜索方法会遇到问题，如果使用SVD从上千篇相似文档中抽取概念，那么同义词会被映射为同一概念。
• 推荐系统。利用SVD从数据中构建一个主题空间，然后再在该空间下计算其相似度。

矩阵分解

• 很多情况下数据中的一小段携带了数据集中的大部分信息。其他信息要么是噪声，要么是毫不相关的信息。矩阵分解可以将原始矩阵表示成新的易于处理的形式，过程类似代数中的因子分解。SVD是最常见的一种矩阵分解技术，SVD将原始数据集Data分解为三个矩阵：U，∑和V^T，这三个矩阵分别是mxm，mxn和nxn。其中矩阵∑只有对角元素，其他元素均为0，并且∑的对角元素是从大到小排列的，这些对角元素称为奇异值，它们对应了原始数据Data的奇异值。这里的奇异值就是矩阵Data*Data^T的特征值的平方根。

• 在某个奇异值的数目（r个）之后，其他奇异值都置为0，这意味着数据集只有r个重要特征，其余特征都是噪声或者冗余特征。

• numpy中的linalg有一个svd方法：U, Sigma, VT = linalg.svd(Data)。注意∑虽然是矩阵，但为了节约空间以array的形式返回。下面的代码(svdRec.py)展示了对样例矩阵求∑的过程。

def loadExData():
        return [[1, 1, 1, 0, 0],
                [2, 2, 2, 0, 0],
                [1, 1, 1, 0, 0],
                [5, 5, 5, 0, 0],
                [1, 1, 0, 2, 2],
                [0, 0, 0, 3, 3],
                [0, 0, 0, 1, 1]]

>>> Data=svdRec.loadExData()
>>> from numpy import *
>>> U, Sigma, VT = linalg.svd(Data)
>>> Sigma
array([  9.72140007e+00,   5.29397912e+00,   6.84226362e-01,
         1.70188300e-16,   5.01684085e-47])

• 从上面返回的∑矩阵可以看出，前三个数值比最后两个值大了很多，因此可以将最后两个值去掉，构成一个3x3的对角矩阵Sig3。若要从U，∑和V^T中构造原始矩阵的近似矩阵，只需要用U的前三列和V^T的前三行。实际操作中，确定要保留的奇异值的个数有两种方法：一是将所有的奇异值map成平方和，之后从前向后叠加，直到累加到总值的90%；二是启发式策略，当矩阵有上万的奇异值时，只保留前3000个，前提是对数据有足够的了解，确保3000个奇异值足够覆盖总平方和的90%。

>>> Sig3 = eye(3)*Sigma[:3]
array([[ 9.72140007,  0.        ,  0.        ],
       [ 0.        ,  5.29397912,  0.        ],
       [ 0.        ,  0.        ,  0.68422636]])
>>> U[:,:3]*mat(Sig3)*VT[:3,:]
matrix([[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
           2.26272993e-16,   2.25622472e-16],
        [  2.00000000e+00,   2.00000000e+00,   2.00000000e+00,
           2.48715978e-16,   2.47198095e-16],
        [  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          -1.18232247e-15,  -1.18297299e-15],
        [  5.00000000e+00,   5.00000000e+00,   5.00000000e+00,
           2.90999863e-16,   2.87530416e-16],
        [  1.00000000e+00,   1.00000000e+00,  -7.77156117e-16,
           2.00000000e+00,   2.00000000e+00],
        [  3.33066907e-16,   6.10622664e-16,  -4.99600361e-16,
           3.00000000e+00,   3.00000000e+00],
        [  1.04083409e-16,   1.87350135e-16,  -1.80411242e-16,
           1.00000000e+00,   1.00000000e+00]])

基于协同过滤的推荐引擎

协同过滤通过将用户和其他用户的数据进行对比来实现推荐。

相似度计算

• 利用不同用户对某一件物品的评分来计算相似度。举例如下表。

  	鳗鱼饭	日式炸鸡	寿司饭	烤牛肉	手撕猪肉
  Jim	2	0	0	4	4
  John	5	5	5	3	3
  Sally	2	4	2	1	2

  可以使用多种方法计算相似度。
  • 欧氏距离：计算手撕牛肉和烤牛肉的距离为sqrt[(4-4)^2+(3-3)^2+(2-1)^2]=1，手撕牛肉和鳗鱼饭的距离为sqrt[(4-2)^2+(3-5)^2+(2-2)^2=2.83。可以使用“相似度=1/(1+距离)”来将相似度控制在0~1之间。
  • 皮尔逊相关系数：该方法较欧氏距离的优势在于，它对用户评分的量级不敏感，例如所有评分都是5分和所有评分都是1分在这里是相同的。numpy中皮尔逊相关系数计算由corrcoef()方法完成，通过0.5+0.5*corrcoef()控制相似度在0~1之间。
  • 余弦相似度：对于两个向量，计算其夹角余弦值来比较相似程度。numpy中提供了linalg.norm()方法用于计算单个向量的2范数（平方和取根）。因为cos在-1~1之间，同样用0.5+0.5*cos来控制相似度在0~1之间。
• 上面几种相似度计算的代码如下，分别为eulidSim，pearsSim和cosSim。

from numpy import *
from numpy import linalg as la

def eulidSim(inA, inB):
    return 1.0/(1.0 + la.norm(inA - inB))

def pearsSim(inA, inB):
    if len(inA) < 3: return 1.0
    return 0.5 + 0.5*corrcoef(inA, inB, rowvar = 0)[0][1]

def cosSim(inA, inB):
    num = float(inA.T * inB)
    denom = la.norm(inA) * la.norm(inB)
    return 0.5+0.5*(num/denom)

>>> myMat = mat(svdRec.loadExData())
>>> svdRec.eulidSim(myMat[:,0], myMat[:,4])
0.13367660240019172
>>> svdRec.cosSim(myMat[:,0], myMat[:,4])
0.54724555912615336
>>> svdRec.pearsSim(myMat[:,0], myMat[:,4])
0.23768619407595826

基于用户还是物品

• 通过计算两种菜肴之间的距离是基于物品的相似度。另一种计算用户距离的方法是基于用户的相似度。在上面的表格示例中，行与行之间的比较是基于用户的相似度，列与列之间的比较是基于物品的相似度。使用哪一种相似度取决于用户或物品的数目。基于X的相似度计算所需的时间会随着X的增长而增长，通常如果用户数目很多并且会不断增长，我们倾向于使用基于物品的相似度计算。

推荐引擎评价

• 通常用于推荐引擎评价的指标是“最小方均根误差（RMSE）”，它计算均方误差的平均值并开根。若评级在1~5分，而RMSE的结果为1.0，说明预测值和用户给出的真实评价差了一分。

餐馆菜肴推荐引擎

给定一个用户，系统会为此用户选择N个最好的推荐菜。整个流程需要做到：寻找用户没有评分的菜肴，在这些没有评分的所有菜肴中，对每种菜计算一个可能的评级分数，即预测用户会对该菜肴做出的评分。最后将评分从高到低排序，返回前N个物品。

• 推荐没有品尝过的菜肴：两个函数，standEst()用于在给定计算相似度方法的前提下，计算用户对某种物品的可能评分；recommend()是推荐引擎，调用standEst函数，并返回前N个最好物品。在standEst的执行过程中，假设要计算用户u对其未打分的菜肴i的可能评分，则需要通过其他物品j和物品i建立联系。扫描所有n个物品，如果用户u对某个物品j有过评分，则寻找所有用户中即对i又对j评分过的用户群体users，根据users们的打分，计算出物品i和物品j的相似度。最后将这个相似度乘以用户u对物品j的评分累加到ratSimTotal变量，将相似度累加到simTotal变量。最后返回ratSimTotal/simTotal就是可能评分。

def standEst(dataMat, user, simMeas, item):
    n = shape(dataMat)[1]
    simTotal = 0.0; ratSimTotal = 0.0
    for j in range(n):
        userRating = dataMat[user,j]
        if userRating == 0: continue
        overLap = nonzero(logical_and(dataMat[:,item].A>0,\
            dataMat[:,j].A>0))[0]
        if len(overLap) == 0: similarity = 0
        else: similarity = simMeas(dataMat[overLap, item],\
            dataMat[overLap, j])
        simTotal += similarity
        ratSimTotal += similarity * userRating
    if simTotal == 0: return 0
    else: return ratSimTotal/simTotal

def recommend(dataMat, user, N=3, simMeas = cosSim, estMethod = standEst):
    unratedItems = nonzero(dataMat[user,:].A==0)[1]
    if len(unratedItems) == 0: return 'you rated everything'
    itemScores = []
    for item in unratedItems:
        estimatedScore = estMethod(dataMat, user, simMeas, item)
        itemScores.append((item, estimatedScore))
    return sorted(itemScores, key=lambda k: k[1], reverse=True)[:N]

>>> myMat = mat(svdRec.loadExData())
>>> myMat[0,1]=myMat[0,0]=myMat[1,0]=myMat[2,0] = 4
>>> myMat[3,3] =2
>>> myMat
matrix([[4, 4, 1, 0, 0],
        [4, 2, 2, 0, 0],
        [4, 1, 1, 0, 0],
        [5, 5, 5, 2, 0],
        [1, 1, 0, 2, 2],
        [0, 0, 0, 3, 3],
        [0, 0, 0, 1, 1]])
>>> svdRec.recommend(myMat,2)
[(4, 2.5), (3, 1.9703483892927431)]
>>> svdRec.recommend(myMat,2,simMeas = svdRec.eulidSim)
[(4, 2.5), (3, 1.9866572968729499)]
>>> svdRec.recommend(myMat,2,simMeas = svdRec.pearsSim)
[(4, 2.5), (3, 2.0)]

• 在代码中加入loadExData2()方法，利用svd处理该矩阵，分析当前矩阵可以发现前三个奇异值就已经占据总能量的90%。

>>> U, Sigma, VT = linalg.svd(mat(svdRec.loadExData2()))
>>> Sigma
array([ 15.77075346,  11.40670395,  11.03044558,   4.84639758,
         3.09292055,   2.58097379,   1.00413543,   0.72817072,
         0.43800353,   0.22082113,   0.07367823])
>>> Sig2 = Sigma**2
>>> sum(Sig2)
541.99999999999966
>>> sum(Sig2)*0.9
487.79999999999973
>>> sum(Sig2[:2])
378.8295595113579
>>> sum(Sig2[:3])
500.50028912757938

• 因为前三个奇异值已经包含了90%的能量，因此可以将11维数据缩减为3维。添加svdEst方法用于简化数据，这里对原书代码做了一点修改，自动提取前90%能量的奇异值。评估相似度时使用的矩阵是一个nxi的矩阵，这里的i是提取的奇异值个数，n是物品数目。不同用户对物品k的评分已经被压缩到第k行的i个数据中，只要计算两个1xi向量的相似度即可。

def svdEst(dataMat, user, simMeas, item):
    n = shape(dataMat)[1]
    simTotal = 0.0; ratSimTotal = 0.0
    U, Sigma, VT = la.svd(dataMat)
    sumTotal = 0
    for singular in Sigma:
        sumTotal += singular**2
    singularMax = 0
    for i in range(len(Sigma)):
        singularMax += Sigma[i]**2
        if singularMax >= sumTotal*0.9:
            break
    SigI = mat(eye(i+1)*Sigma[:i+1])
    xformedItems = dataMat.T * U[:,:i+1] * SigI.I
        # n x m     m x i    i x i   =>     n x i
    for j in range(n):
        userRating = dataMat[user,j]
        if userRating == 0 or j == item: continue
        similarity = simMeas(xformedItems[item,:].T,\
            xformedItems[j,:].T)
        print 'the %d and %d similarity is: %f' % (item, j, similarity)
        simTotal += similarity
        ratSimTotal += similarity * userRating
    if simTotal == 0: return 0
    else: return ratSimTotal/simTotal

>>> myMat = mat(svdRec.loadExData2())
>>> svdRec.recommend(myMat, 1, N=3, estMethod = svdRec.svdEst)
the 0 and 3 similarity is: 0.485722
the 0 and 5 similarity is: 0.486944
....
[(6, 3.3329499901459845), (9, 3.3315447178728395), (4, 3.3314474877128624)]
>>> svdRec.recommend(myMat, 1, simMeas = svdRec.pearsSim)
[(6, 3.3333333333333335), (9, 3.3333333333333335), (0, 3.0)]

• SVD在大数据集上会显著降低程序速度，可以仅在程序调入运行时离线加载一次。其中计算出的相似度是物品和物品之间的相似度，不同的用户也可以重复使用，因此可以将计算结果离线。当出现冷启动问题时，可以将推荐看成搜索问题，为物品添加标签，使用基于内容的推荐。

基于SVD的图像压缩

• 用一个32*32（1024像素）的矩阵表示一个图像，通过svd对该矩阵降维，实现压缩重构。运行结果太长，但基本还原了原来的矩阵（可以观察出部分像素点不同）。U和V^T都是32x2的矩阵，有两个奇异值，因此总数字数目为64x2+2=130，获得了几乎10倍的压缩比。

def printMat(inMat, thresh=0.8):
    for i in range(32):
        for k in range(32):
            if float(inMat[i,k]) > thresh:
                print 1,
            else: print 0,
        print ''

def imgCompress(numSV=3, thresh=0.8):
    myl = []
    for line in open('0_5.txt').readlines():
        newRow = []
        for i in range(32):
            newRow.append(int(line[i]))
        myl.append(newRow)
    myMat = mat(myl)
    print "****original matrix******"
    printMat(myMat, thresh)
    U,Sigma,VT = la.svd(myMat)
    SigRecon = mat(zeros((numSV, numSV)))
    for k in range(numSV):#construct diagonal matrix from vector
        SigRecon[k,k] = Sigma[k]
    reconMat = U[:,:numSV]*SigRecon*VT[:numSV,:]
    print "****reconstructed matrix using %d singular values******" % numSV
    printMat(reconMat, thresh)

SVD总结

SVD可以逼近矩阵并从中提取重要特征。通过保留矩阵80%~90%的能量，可以去除噪声，保留重要特征，在稀疏矩阵的压缩上作用显著。SVD的其中一个应用是推荐引擎，协同过滤基于用户喜好或者行为数据来推荐，核心是相似度计算方法，SVD可以将高维用户群体降维成少数奇异值，提高推荐引擎效果。在大规模数据集上，SVD的冗余计算会占用过多时间，可以通过离线方式SVD分解和相似度计算。

---

参考文献： 《机器学习实战 - 美Peter Harrington》

原创作品，允许转载，转载时无需告知，但请务必以超链接形式标明文章原始出处(https://forec.github.io/2016/02/26/machinelearning14/) 、作者信息（Forec）和本声明。