《机器学习实战》之十四——利用SVD简化数据

一、前言

  奇异值分解是一个有着很明显的物理意义的一种方法，它可以将一个比较复杂的矩阵用更小更简单的几个子矩阵的相乘来表示，这些小矩阵描述的是矩阵的重要的特性。就像是描述一个人一样，给别人描述说这个人长得浓眉大眼，方脸，络腮胡，而且带个黑框的眼镜，这样寥寥的几个特征，就让别人脑海里面就有一个较为清楚的认识，实际上，人脸上的特征是有着无数种的，之所以能这么描述，是因为人天生就有着非常好的抽取重要特征的能力，让机器学会抽取重要的特征，SVD是一个重要的方法。

  通过SVD对数据的处理，我们可以使用小得多的数据集来表示原始数据集，这样  做实际上是去除了噪声和冗余信息，以此达到了优化数据、提高结果的目的。

  在机器学习领域，有相当多的应用与奇异值都可以扯上关系，比如做 feature reduction 的PCA，做数据压缩（以图像压缩为代表）的算法，还有做搜索引擎语义层次检索的LSI（Latent Semantic Indexing）等。

二、SVD的应用

1、隐性语义索引

  最早的SVD应用之一就是信息检索，我们称利用SVD的方法为隐性语义索引（Latent Semantic Indexing，LSI）或 隐形语义分析（Latent Semantic Analysis，LSA）。

  在LSI中，一个矩阵是由 文档和词语组成的。当我们在该矩阵上应用SVD时，就会构建出多个奇异值。这些奇异值代表了文档的 概念或主题，这一特点可以用于更高效的文档搜索。在词语拼写错误时，只基于词语存在与否的简单搜索方法会遇到问题。简单搜索的另一个问题就是同义词的使用。这就是说，当我们查找一个词时，其同义词所在的文档可能并不会匹配上。如果我们从上千篇相似的文档中抽取出概念，那么同义词就会映射为同一个概念。

2、推荐系统

  SVD的另一个应用就是推荐系统，较为先进的推荐系统先利用SVD从数据中构建一个主题空间，然后再在该空间下计算相似度，以此提高推荐的效果。

三、SVD的原理

  特征值分解和奇异值分解两者有着很紧密的关系，我在接下来会谈到，特征值分解和奇异值分解的目的都是一样，就是提取出一个矩阵最重要的特征。先谈谈特征值分解吧：

1、特征值分解

  如果说一个向量 $v$ 是方阵A的特征向量，将一定可以表示成下面的形式：
$$Av=\lambda v$$

  这时候 $\lambda$ 就被称为特征向量 $v$ 对应的特征值，一个矩阵的一组特征向量是一组正交向量。特征值分解是将一个矩阵分解成下面的形式：

$$A=Q\Sigma Q^1$$

  其中Q是这个矩阵A的特征向量组成的矩阵，$\Sigma$ 是一个对角阵，每一个对角线上的元素就是一个特征值。首先，要明确的是，一个矩阵其实就是一个线性变换，因为一个矩阵乘以一个向量后得到的向量，其实就相当于将这个向量进行了线性变换。比如说下面的一个矩阵：

$$M=\left[\begin{array}{cc}3& 0\\ 0& 1\end{array}\right]$$
它其实对应的线性变换是下面的形式：

因为这个矩阵M乘以一个向量(x,y)的结果是：
$$\left[\begin{array}{cc}3& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}3x\\ y\end{array}\right]$$
  上面的矩阵是对称的，所以这个变换是一个对x，y轴的方向一个拉伸变换（每一个对角线上的元素将会对一个维度进行拉伸变换，当值>1时，是拉长，当值<1时时缩短），当矩阵不是对称的时候，假如说矩阵是下面的样子：
$$M=\left[\begin{array}{cc}3& 0\\ 0& 1\end{array}\right]$$

它所描述的变换是下面的样子：

  这其实是在平面上对一个轴进行的拉伸变换（如蓝色的箭头所示），在图中，蓝色的箭头是一个最主要的变化方向（变化方向可能有不止一个），如果我们想要描述好一个变换，那我们就描述好这个变换主要的变化方向就好了。反过头来看看之前特征值分解的式子，分解得到的 $\Sigma$ 矩阵是一个对角阵，里面的特征值是由大到小排列的，这些特征值所对应的特征向量就是描述这个矩阵变化方向（从主要的变化到次要的变化排列）

  当矩阵是高维的情况下，那么这个矩阵就是高维空间下的一个线性变换，这个线性变化可能没法通过图片来表示，但是可以想象，这个变换也同样有很多的变换方向，我们通过特征值分解得到的前N个特征向量，那么就对应了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向，就可以近似这个矩阵（变换）。也就是之前说的：提取这个矩阵最重要的特征。

  总结一下，特征值分解可以得到特征值与特征向量，特征值表示的是这个特征到底有多重要，而特征向量表示这个特征是什么，可以将每一个特征向量理解为一个线性的子空间，我们可以利用这些线性的子空间干很多的事情。不过，特征值分解也有很多的局限，比如说变换的矩阵必须是方阵。

2、奇异值分解

  特征值分解是一个提取矩阵特征很不错的方法，但是它只是对 方阵 而言的，在现实的世界中，我们看到的大部分矩阵都不是方阵，比如说有M个学生，每个学生有N科成绩，这样形成的一个M * N的矩阵就不可能是方阵，我们怎样才能描述这样普通的矩阵的重要特征呢？奇异值分解可以用来干这个事情，奇异值分解是一个能适用于任意的矩阵的一种分解的方法：
$$A=U\Sigma V^T$$
  假设A是一个M * N的矩阵，那么得到的U是一个M * M的方阵（里面的向量是正交的，U里面的向量称为左奇异向量），$\Sigma$ 是一个M * N的矩阵（除了对角线的元素都是0，对角线上的元素称为奇异值），$V^T$是一个N * N的矩阵，里面的向量也是正交的，V里面的向量称为右奇异向量），从图片来反映几个相乘的矩阵的大小可得下面的图片

  那么奇异值和特征值是怎么对应起来的呢？首先，我们将一个矩阵 $A^T\ast A$，将会得到一个方阵，我们用这个方阵求特征值可以得到：

$$(A^{T}A)v_i={\lambda }_i{v}_i$$

  这里得到的$v$，就是我们上面的右奇异向量。此外我们还可以得到：

$$\begin{gathered} {\sigma }_i=\sqrt{{\lambda }_i} \\ {u}_i=\frac{1}{{\sigma }_i}A{v}_i \end{gathered}$$

  这里的 $\sigma$ 就是上面说的奇异值，$u$ 就是上面说的左奇异向量。奇异值$\sigma$跟特征值类似，在矩阵 $\Sigma$ 中也是从大到小排列，而且 $\sigma$ 的减少特别的快，在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上了。也就是说，我们也可以用前 $r$ 大的奇异值来近似描述矩阵，这里定义一下部分奇异值分解： 即前$r$个非零奇异值对应的奇异向量代表了矩阵的主要特征,

$$A_{m\times n}\approx U_{m\times r}{\Sigma }_{r\times r}{V}_{r\times n}^T$$

  $r$ 是一个远小于m、n的数，这样矩阵的乘法看起来像是下面的样子：

   右边的三个矩阵相乘的结果将会是一个接近于A的矩阵，在这儿，$r$ 越接近于n，则相乘的结果越接近于A。而这三个矩阵的面积之和（在存储观点来说，矩阵面积越小，存储量就越小）要远远小于原始的矩阵A，我们如果想要压缩空间来表示原矩阵A，我们存下这里的三个矩阵：$U$、$\Sigma$、$V$就好了。

四、利用python实现SVD

  如果SVD确实那么好，那么如何实现它呢？SVD实现了相关额线性代数。NumPy有一个称为linalg的线性代数工具箱，接下来，我们了解一下如何利用该工具箱实现如下矩阵的SVD处理。

$$\left[\begin{array}{cc}1& 1\\ 1& 7\end{array}\right]$$
新建一个svdRec.py文件，在其中写入如下的python代码：

# -*- coding: utf-8 -*-

import numpy as np

data = [[1,1],[1,7]]
U,Sigma,VT = np.linalg.svd(np.mat(data))
print("U:\n",U)
print("Sigma:\n",Sigma)
print("VT:\n",VT)

运行之后结果如下：

我们从结果中注意到，矩阵Sigma是以行向量array([7.16227766, 0.83772234]) 返回，而非如下矩阵
$$\left[\begin{array}{cc}7.16227766& 0.0\\ 0.0& 0.83772234\end{array}\right]$$

由于矩阵除了对角元素其他均为0，因此这种仅返回对角元素的方式能够节省空间，这就是由Numpy内部的机制产生的。我们所要记住的是，一旦看到Sigma就要知道它是一个矩阵。

打开svdRec.py文件，在其中写入如下的python代码：

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]]
    
if __name__ == "__main__":
    Data = loadExData()
    U,Sigma,VT = np.linalg.svd(Data)
    print("Sigma:\n",Sigma)

运行结果如下：

  这里要注意的是在不同的机器上产生的结果可能会稍有不同，尤其是最后两个值，因为他们太小了。而前3个数值结果应该是一样的。

  接下来，我们的原始数据集就可以用如下的结果来近似。
$$Dat{a}_{m\times n}\approx U_{m\times 3}{\Sigma }_{3\times 3}{V}_{3\times n}^T$$
  我们用Sigma的前3个值来重构原始矩阵的近似矩阵。图解如下：

代码如下：

import numpy as np

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]]
    
if __name__ == "__main__":
    Data = loadExData()
    U,Sigma,VT = np.linalg.svd(Data)
    print("Sigma:\n",Sigma)
    Sig3 = np.mat([[Sigma[0],0,0],[0,Sigma[1],0],[0,0,Sigma[2]]]) #构建3×3的Sigma矩阵
    Data_ = U[:,:3] * Sig3 * VT[:3,:]
#    print("近似矩阵：\n",Data_.round(2))
    Data_.round(2)

运行可得到近似矩阵Data_的结果如下：

  现在我们通过三个矩阵$U_{m\times 3}$、${\Sigma }_{3\times 3}$、$V_{3\times n}^T$ 对原始矩阵进行了近似。从结果来看，与原始矩阵是一样的，可见，我们可以用一个小很多的矩阵来表示一个大矩阵。

1、基于协同过滤的推荐引擎

  近年来，推荐引擎对因特网用户而言已经不是什么新鲜事物了。有很多方法可以实现推荐功能，这里我们使用的是一种称为协同过滤（collaborative filtering）的方法。
  大致有两种方法来实现：

• 基于用户 (user-based) 的协作型过滤
• 基于物品 (item-based) 的协作型过滤

  两种方法大致相同，但是在不同的环境下，使用最佳的方法能最大化的提升算法的效果。如下图（后面的示例数据）所示，对两样商品直接的距离进行计算，这称为基于物品的相似度。而对行与行（用户之间）进行距离的计算，这称为基于用户的相似度。到底该选用那种方法呢？这取决与用户或物品的数量，基于物品相似度的计算时间会随着物品数量的增加而增加。基于用户相似度则取决于用户数量，例如：一个最大的商店拥有大概100000种商品，而它的用户可能有500000人，这时选择基于物品相似度可能效果好很多。

2、相似度计算方法：

相似度的计算方法常用的有三种方法：

• 欧氏距离

• 皮尔逊相关系数（Pearson correlation）

• 余弦相似度（cosine similarity）
  计算余弦相似度，我们采用的两个向量A和B夹角的余弦相似度。如果夹角为90度，则相似度为0；如果两个向量的方向相同，则相似度为1.0。其定义如下：
  $$cos\theta =\frac{A\cdot B}{\Vert A\Vert \Vert B\Vert }$$
  其中，$\Vert A\Vert$、$\Vert B\Vert$表示向量A、B的2范数。向量 [4,2,2] 的2范数为
  $$\sqrt{4^2+3^2+2^2}$$

接下来，我们打开svdRec.py文件，并在其中写入相似度计算的代码：

def ecludSim(inA, inB):
    return 1.0/(1.0 + np.linalg.norm(inA - inB)) #返回结果已转换到0，1   1最大相似，0最小相似

def pearsSim(inA, inB):
    if len(inA) <3: 
        return 1.0
    return 0.5+0.5*np.corrcoef(inA, inB, rowvar=0)[0][1] # 使用0.5+0.5*x 将-1，1 转为 0，1

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

五、示例：餐馆菜肴推荐引擎

  现在我们就开始构建一个推荐引擎，该推荐引擎关注的是餐馆食物的推荐。假设一个人在家决定外出吃饭，但是他并不知道该到哪儿去吃饭，该点什么菜。我们这个推荐系统可以帮助他做到这两点。

1、推荐未尝过的菜肴

  推荐系统的工作过程是：给定一个用户，系统会为此用户返回N个最好的推荐菜。为了实现这一点，则需要我们做到：

• （1）寻找用户没有评级过的菜肴，即在用户-物品矩阵中的0值
• （2）在用户没有评级的所有物品总，对每个物品预计一个可能的评级分数。这就是说，我们认为如果用户对该物品打分的话能打多少分（这就是相似度计算的初衷）
• （3）对这些物品的评分从高到低进行排序，返回前N个物品

  接下来，我们打开svdRec.py文件，并在其中写入如下的代码：

"""
函数说明 ：用户对物品的估计的评分函数
Parameters：
    dataMat: 数据集
    user: 用户编号
    simMeas: 相似度度量方法
    item: 物品编号
Returns：
    用户对物品的估计评分值
"""
def standEst(dataMat, user, simMeas, item):
    n = np.shape(dataMat)[1] #计算物品的数目
    simTotal = 0.0 
    ratSimTotal = 0.0
    for j in range(n): #循环遍历所有物品
        userRating = dataMat[user, j]  #获取用户对每个物品j的评分
        if userRating == 0: #如果用户对物品没有评分则跳过该物品
            continue
        ## 获得相比较的两个物品中不为0的重合元素的行号
        overLap = np.nonzero(np.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])
       # 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 # 归一化评分 使其处于0-5（评级）之间
"""
函数说明 ：物品推荐函数
Parameters：
    dataMat:数据集
    user: 用户编号
    N: 推荐数目
    simMeas: 相似度度量方法
    estMethod: 估计得分的方法
Returns：
    包含（推荐的物品号，得分）的列表，且该列表按照得分进行降序排列
"""
def recommend(dataMat, user, N=3, simMeas=cosSim, estMethod=standEst):
    unratedItems = np.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 jj:jj[1], reverse=True)[:N] #按照得分进行降序排列


if __name__ == "__main__":
#    Data = loadExData()
#    U,Sigma,VT = np.linalg.svd(Data)
#    print("Sigma:\n",Sigma)
#    Sig3 = np.mat([[Sigma[0],0,0],[0,Sigma[1],0],[0,0,Sigma[2]]])
#    Data_ = U[:,:3] * Sig3 * VT[:3,:]
##    print("近似矩阵：\n",Data_.round(2))
#    Data_.round(2)

    data = [[4,4,0,2,2],
            [4,0,0,3,3],
            [4,0,0,1,1],
            [1,1,1,2,0],
            [2,2,2,0,0],
            [1,1,1,0,0],
            [5,5,5,0,0]]
    myMat = np.mat(data)
    result1 = recommend(myMat, 2)
    print("余弦相似计算结果：\n",result1)
    result2 = recommend(myMat,2, simMeas=ecludSim)
    print("欧氏距离计算结果：\n",result2)
    result3 = recommend(myMat, 2, simMeas=pearsSim)
    print("皮尔逊相似度计算结果：\n",result3)
    

运行结果如下图所示：

2、利用SVD提高推荐的效果

  实际的数据集会比我们用于展示的recommend（）函数功能的myMat矩阵 稀疏 的多。下面我们计算一个比较稀疏的矩阵的SVD来看看到底需要多少维特征。

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

if __name__ == "__main__":
    myMat = np.mat(loadExData2())
    U, Sigma, VT = np.linalg.svd(myMat)
    print("The length of Sigma:\n", len(Sigma))

  运行结果，可以的Sigma的长度是11，接下来，我们看一下Sigma中的前几个元素可以达到总值的90%。打开svdRec.py文件，在main函数中继续写入如下的代码：

	Sig2 = Sigma**2  #对Sigma中的值求平方
    Energy = sum(Sig2)
    print("总能量：", Energy)
    print("90%的总能量：",Energy*0.9)
    print("前两个元素的能量：", sum(Sig2[:2]))
    print("前三个元素的能量：", sum(Sig2[:3]))

运行结果如下：

  从结果我们可以看出，前三个元素的能量已经达到了总能量的90%了。其实，我们可以用一个3维的矩阵来代替11维的Sigma矩阵。

  这里要注意的是，总能量的计算是通过对Sigma中的值求平方后再求和得到的。为什么要求平方呢？我想这里是怕Sigma中的值有负值，那样总能量计算就没有可比性了。

  接下来，我们利用SVD将所有的菜肴映射到一个低维的空间中去。在低维空间下，可以利用前面相同的相似度计算方法来进行推荐。打开svdRec.py文件，继续写入如下的代码：

"""
函数说明 ：基于SVD降维后用户对物品估计评分函数
Parameters：
    dataMat: 数据集
    user: 用户编号
    simMeas: 相似度度量方法
    item: 物品编号
Returns：
    用户对物品的估计评分值
"""
def svdEst(dataMat, user, simMeas, item):
    n = np.shape(dataMat)[1]
    simTotal = 0.0
    ratSimTotal = 0.0
    U,Sigma, VT = np.linalg.svd(dataMat)
    Sig4 = np.mat(np.eye(4)*Sigma[:4]) #将4维的奇异值向量转换为奇异值对角矩阵
    
    # 降维方法 通过U矩阵将物品转换到低维空间中 （商品数行x选用奇异值列）
    xformedItems = dataMat.T *U[:,:4] * Sig4.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

if __name__ == "__main__":
    myMat = np.mat(loadExData2())
    recom1 = recommend(myMat, 1, estMethod=svdEst)
    print("推荐结果1：\n",recom1)
    recom2 = recommend(myMat, 1, estMethod=svdEst, simMeas=pearsSim)
    print("推荐结果2：\n",recom2)

运行结果如下图所示：


  在函数svdEst(）中之所以使用Sig4（将原始数据降到4维），是因为通过预先计算得到能满足90%的奇异值能量的是前4个奇异值。

  在函数svdEst() 中使用SVD方法，将数据集映射到低维的空间中，再做运算。其中的 xformedItems = dataMat.T*U[:,:4]*Sig4.I 可能不是很好理解，它就是SVD的降维步骤，通过U矩阵和Sig4逆矩阵将商品转换到低维空间（得到 商品行，选用奇异值列）。

六、示例：基于SVD的图像压缩

  接下来，我们要做的是一个很好的关于如何将SVD应用于图像压缩的例子。通过可视化的方式，该例子使得我们很容易就能看到SVD对数据近似的效果。

  这里不采用书中的例子来讲解，因为无趣所以这里换作我们的真实的照片来做一个简单的SVD图片压缩作为一个示例：

首先放上男神图片：

  基于SVD图片压缩原理其实很简单，图片其实就是数字矩阵，通过SVD将该矩阵降维，只使用其中的重要特征来表示该图片从而达到了压缩的目的。

# -*- coding: utf-8 -*-

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

"""
函数说明：svd图片压缩函数
Parameters：
    k：前k个
    pic_array:图片数据矩阵
Returns：
    new_pic: 压缩后的图片数据
    size: 压缩后的图片大小
"""
def pic_compress(k, pic_array):
    U, Sigma, VT = np.linalg.svd(pic_array) #用svd分解图片数据矩阵
    SigK = np.eye(k) * Sigma[: k] # k维的sigma对角矩阵
    new_pic = U[:, :k] * SigK * VT[:k, :]  # 重构的图片矩阵
    size = U.shape[0] * k + SigK.shape[0] * SigK.shape[1] + k * VT.shape[1]  # 压缩后大小
    return new_pic, size



ori_img = np.array(ndimage.imread("wed.jpg", flatten=True))  #读入原始图片
new_img, size = pic_compress(30, ori_img) #取前30个奇异值来还原原图像
print("original size:" +str(ori_img.shape[0])+"*"+str(ori_img.shape[1])+"="+str(ori_img.shape[0] * ori_img.shape[1]) )
print("compress size:" + str(size))
fig, ax = plt.subplots(1, 2)
ax[0].imshow(ori_img,cmap='gray')
ax[0].set_title("before compress")
ax[1].imshow(new_img,cmap='gray')
ax[1].set_title("after compress")
plt.show()

当取前10个奇异值，运行结果如下所示：

当取前30个奇异值，运行结果如下所示：

当取前50个奇异值，运行结果如下所示：

  比较压缩前后两张图片，取前50个奇异值时，图片已经看着很清晰了。接下来，我们看一下压缩前后图片的大小。

  原图片为870x870，保存像素点值为870x870 = 756900，使用SVD，取前50个奇异值则变为：89500
存储量大大减小，仅50个奇异值就已经能很好的反应原数据了。

  为什么选取前50个奇异值，其实我们可以根据前面介绍的方法进行预先的测试。

  值得一提的是，奇异值从大到小衰减得特别快，在很多情况下，前 10% 甚至 1% 的奇异值的和就占了全部的奇异值之和的 99% 以上了。这对于数据压缩来说是个好事。下面这张图展示了本例中奇异值和奇异值累加的分布。

上述绘制图的代码如下：

def PlotPic(pic_array):
    U, Sigma, VT = np.linalg.svd(pic_array) #用svd分解图片数据矩阵
    cnt = []
    temp = 0
    m = len(Sigma)
    for i in range(m):
        temp += Sigma[i]
        cnt.append(temp)
    plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
    plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
    
    plt.plot(Sigma,'r')  
    plt.plot(cnt,'b')
    label = ["Sigular Value", "Sum of Sigular Value"]  
    plt.legend(label, loc = 5, ncol = 1)
    plt.xticks(np.linspace(0,1000,6))  #设置x坐标轴的刻度，范围0-20，生成共5个点
    plt.show()

ori_img = np.array(ndimage.imread("wed.jpg", flatten=True))  #读入原始图片
PlotPic(ori_img)

SVD两个最重要的计算步骤下：

• 数据集降维： 这里的sigma为对角矩阵（需要利用原来svd返回的sigma向量构建矩阵，构建需要使用count这个值）。U为svd返回的左奇异矩阵，count为我们指定的多少个奇异值，这也是sigma矩阵的维数。
• 重构数据集： 这里的sigma同样为对角矩阵（需要利用原来svd返回的sigma向量构建矩阵，构建需要使用count这个值），VT为svd返回的右奇异矩阵，count为我们指定的多少个奇异值（可以按能量90%规则选取）。

七、总结

• SVD是一种强大的降维工具，我们可以用SVD来逼近矩阵并从中提取重要特征。去除噪声
• 协同过滤是一种基于用户喜好或行为数据的推荐的实现方法。通过在低维空间下计算相似度，SVD提高了推荐系统引擎的效果。
• 在大规模数据集上，SVD的计算和推荐可能是一个很困难的工程问题。通过离线方式来进行SVD分解和相似度计算，是一种减少冗余计算和推荐所需时间的办法。

参考资料

【1】奇异值（SVD）原理详解：https://blog.csdn.net/xiaocong1990/article/details/54909126
【2】《机器学习实战》第14章 利用SVD简化数据
【3】https://www.zhihu.com/question/22237507/answer/53804902
【4】https://blog.csdn.net/wang454592297/article/details/80999644
【5】如何让奇异值分解(SVD)变得不“奇异”？ ： http://redstonewill.com/1529/