奇异值分解实验

  • 奇异值分解
  • 低秩近似
  • 工程应用:图像压缩
  • 工程应用:推荐系统



 


奇异值分解

只有方阵(行数等于列数)才能做特征值分解,非方阵可不可以分解为 奇异值分解实验:图像压缩与推荐系统_奇异值分解

这种方式是【奇异值分解】,这种方法大学里并不学。

因为本科的线性代数主要研究方阵(除了线性系统),所以大学里并没有介绍非方阵的奇异值分解(奇异值分解实验:图像压缩与推荐系统_特征值_02),奇异值分解在数据降维、语义分析、图像等领域都有十分广泛的应用,比如 奇异值分解实验:图像压缩与推荐系统_工程应用_03

举个荔枝,演示一下非方阵的分解步骤。

  • 奇异值分解实验:图像压缩与推荐系统_奇异值分解_04
  1. 求出矩阵 奇异值分解实验:图像压缩与推荐系统_特征值_05奇异值分解实验:图像压缩与推荐系统_工程应用_06
    奇异值分解实验:图像压缩与推荐系统_奇异值分解_07
     
    奇异值分解实验:图像压缩与推荐系统_奇异值分解_08
    发现 奇异值分解实验:图像压缩与推荐系统_特征值_05奇异值分解实验:图像压缩与推荐系统_工程应用_06
     
  2. 因此,分别求出 奇异值分解实验:图像压缩与推荐系统_特征值_05奇异值分解实验:图像压缩与推荐系统_工程应用_06
    奇异值分解实验:图像压缩与推荐系统_特征值_05 有俩组,奇异值分解实验:图像压缩与推荐系统_奇异值分解_14
     
    奇异值分解实验:图像压缩与推荐系统_工程应用_06 有三组,奇异值分解实验:图像压缩与推荐系统_奇异值分解_16
     
  3. 奇异值分解实验:图像压缩与推荐系统_工程应用_06 的特征向量横向拼成矩阵 奇异值分解实验:图像压缩与推荐系统_奇异值分解_18
    奇异值分解实验:图像压缩与推荐系统_特征值_19
    奇异值分解实验:图像压缩与推荐系统_奇异值分解_18
     
  4. 奇异值分解实验:图像压缩与推荐系统_特征值_05奇异值分解实验:图像压缩与推荐系统_工程应用_06 的相同特征值开方,拼成矩阵奇异值分解实验:图像压缩与推荐系统_特征值_23
     
  5. 奇异值分解实验:图像压缩与推荐系统_工程应用_24

  6. 拼成矩阵奇异值分解实验:图像压缩与推荐系统_特征值_23奇异值分解实验:图像压缩与推荐系统_奇异值分解_26
    P.S. 奇异值分解实验:图像压缩与推荐系统_奇异值分解_27
     
  7. 奇异值分解实验:图像压缩与推荐系统_特征值_05 的特征向量横向拼成一个矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_29
    奇异值分解实验:图像压缩与推荐系统_工程应用_30
    P.S. 奇异值分解实验:图像压缩与推荐系统_奇异值分解_31 也是正交矩阵,因为TA的列向量俩俩正交,且是单位向量。
     
  8. 最后,将 奇异值分解实验:图像压缩与推荐系统_特征值_32 发现结果等于原矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_33,说明 矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_33 可以分解为奇异值分解实验:图像压缩与推荐系统_特征值_35
    奇异值分解实验:图像压缩与推荐系统_特征值_36

虽然这只是一个 奇异值分解实验:图像压缩与推荐系统_奇异值分解_37 的矩阵的分解过程,但推广到 奇异值分解实验:图像压缩与推荐系统_工程应用_38

  • 奇异值分解:奇异值分解实验:图像压缩与推荐系统_工程应用_39
    奇异值分解实验:图像压缩与推荐系统_奇异值分解_40, 矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_41 的尺寸是 奇异值分解实验:图像压缩与推荐系统_特征值_42
    奇异值分解实验:图像压缩与推荐系统_特征值_43,矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_44 的尺寸是 奇异值分解实验:图像压缩与推荐系统_特征值_45,其列向量称为 矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_41
    奇异值分解实验:图像压缩与推荐系统_特征值_47,矩阵奇异值分解实验:图像压缩与推荐系统_特征值_48 的尺寸是 奇异值分解实验:图像压缩与推荐系统_特征值_42,其对角线上的值称为 矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_41
    奇异值分解实验:图像压缩与推荐系统_奇异值分解_51, 矩阵奇异值分解实验:图像压缩与推荐系统_工程应用_52 的尺寸是 奇异值分解实验:图像压缩与推荐系统_特征值_53,其列向量称为 矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_41

奇异值分解实验:图像压缩与推荐系统_工程应用_55

 


低秩近似

我们可以把矩阵看成一种变换,把矩阵乘法当成线性变换时,找出变换矩阵的特征值和特征向量,实际上就是找出变换矩阵的主要变换方向。

也可以说是,特征值和特征向量代表了一个方阵的【固有信息】。

特征值分解是奇异值分解的特例,特征值分解只能分解方阵,奇异值分解可以分解任意形状的矩阵。

因此,奇异值及奇异向量可以说是代表了一个 奇异值分解实验:图像压缩与推荐系统_工程应用_38

奇异值越大,代表的信息就越多。

另外,如果我们在奇异值矩阵 奇异值分解实验:图像压缩与推荐系统_特征值_57

很多情况下,前 奇异值分解实验:图像压缩与推荐系统_特征值_58 个奇异值的和就占了全部奇异值之和的 奇异值分解实验:图像压缩与推荐系统_奇异值分解_59%。

也就是说,前 奇异值分解实验:图像压缩与推荐系统_特征值_58

所以,可以用 最大的 奇异值分解实验:图像压缩与推荐系统_特征值_58

这种理论,就被称为【低秩近似】。

来看一个 奇异值分解实验:图像压缩与推荐系统_工程应用_62

  • 奇异值分解实验:图像压缩与推荐系统_奇异值分解_63

    奇异值分解:奇异值分解实验:图像压缩与推荐系统_工程应用_39,看中间的奇异值分解实验:图像压缩与推荐系统_特征值_48矩阵: 奇异值分解实验:图像压缩与推荐系统_奇异值分解_66

    发现奇异值有 奇异值分解实验:图像压缩与推荐系统_工程应用_67 个,分别是 奇异值分解实验:图像压缩与推荐系统_工程应用_68,我们只用前面俩个奇异值来算一下,占总奇异值的比例:
  • 奇异值分解实验:图像压缩与推荐系统_特征值_69%

 
这个比例很大了,所以我们可以认为前面俩个奇异值,及其对应的左右奇异向量足以代表原来的矩阵。

奇异值分解实验:图像压缩与推荐系统_特征值_70

把这些部分截取下来:

  • 奇异值分解实验:图像压缩与推荐系统_奇异值分解_71

将截取下来的矩阵相乘,就低秩近似原来的数据矩阵,对比俩个矩阵可以发现,俩者数值非常接近。

  • 原始矩阵: 奇异值分解实验:图像压缩与推荐系统_工程应用_72 近似矩阵:奇异值分解实验:图像压缩与推荐系统_奇异值分解_73

 


工程应用:图像压缩

基于奇异值分解的【低秩近似】理论在工程中有广泛的应用,比如图像压缩。

图像本来就是一个矩阵,比如下面的图片:

奇异值分解实验:图像压缩与推荐系统_特征值_74


假设这张图片的尺寸是 奇异值分解实验:图像压缩与推荐系统_工程应用_75,就需要 奇异值分解实验:图像压缩与推荐系统_工程应用_75

就这样一张不大的灰度图,都要将近 奇异值分解实验:图像压缩与推荐系统_奇异值分解_77 万字节(奇异值分解实验:图像压缩与推荐系统_特征值_78)存储,要是某个应用是图片为主的,那您可以想象应用会有多大。

图像压缩主要有俩个好处:

  • 存储空间会小很多
  • 方便网络传输

我们可以用 奇异值分解 来压缩图像,算法就是【低秩近似】理论。

图像压缩有俩部分:

  • 压缩图像
    奇异值分解实验:图像压缩与推荐系统_特征值_42 图像矩阵做奇异值分解,得到 奇异值分解实验:图像压缩与推荐系统_特征值_80

    选取前 奇异值分解实验:图像压缩与推荐系统_工程应用_81 大的奇异值(奇异值分解实验:图像压缩与推荐系统_特征值_82),按照【低秩近似】理论,对 奇异值分解实验:图像压缩与推荐系统_特征值_80 做截取,得到 奇异值分解实验:图像压缩与推荐系统_奇异值分解_84

奇异值分解实验:图像压缩与推荐系统_奇异值分解_85

存储或者传输新的 奇异值分解实验:图像压缩与推荐系统_奇异值分解_84,新的矩阵不一定比原来的小,一定要选取一个恰当的 奇异值分解实验:图像压缩与推荐系统_工程应用_81

  • 图像重构
    奇异值分解实验:图像压缩与推荐系统_奇异值分解_84

完整代码:

import cv2
# cv 库用来读取图片
import numpy as np
import matplotlib.pyplot as plt
# plt 库显示图片

''' @para: c 是保留奇异值(奇异向量)个数占总个数比例 '''
def imgCompress(c, img):
    # 1. 图像压缩(SVD), 返回的是分解后的 3 个矩阵
    U, sigma, VT = np.linalg.svd(img)
    k = int(c * img.shape[1])                                   # .shape[1] 是列数,也是奇异值的个数
    sig = np.eye(k) * sigma[ :k]                                # sigma矩阵截取,构造新的奇异值矩阵,也是一个对角矩阵

    # 2. 图像重构
    res_img = (U[:, :k] * sig) * VT[:k, :]                      # U、VT矩阵截取,并相乘
    size = U.shape[0] * k + sig.shape[0] * sig.shape[1] + k * VT.shape[1]
    # 压缩后的数据量 = 截取后的(U的大小 + sigma的大小 + VT的大小)

    return res_img, size;


if __name__ == '__main__':
    # 1.读取待压缩的图像
    img_path = input("图片路径:>  ")
    ori_img = np.mat(cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)) #  cv2.IMREAD_GRAYSCALE:以灰度图方式读取

    # 2.图像压缩(含重构)
    res_img, size = imgCompress(0.1, ori_img)                    # 0.1 只保留前 10% 的奇异值,比例越小,压缩的越厉害,重构的图片就越模糊

    print("压缩前图像大小:>  ", str(ori_img.shape[0] * ori_img.shape[1]))
    print("压缩后图像大小:>  ", str(size))

    # 显示图像(对比)
    fig, ax = plt.subplots(1, 2)
    ax[0].imshow(ori_img, cmap='gray')
    ax[0].set_title("before compress")
    ax[1].imshow(res_img, cmap='gray')
    ax[1].set_title("after compress")
    plt.show()

运行结果:

  • 压缩前图像大小:> 50292
  • 压缩后图像大小:> 8949哈哈,整整压缩了一个量级(奇异值分解实验:图像压缩与推荐系统_奇异值分解_89)。

显示图片:

奇异值分解实验:图像压缩与推荐系统_特征值_90

 


工程应用:推荐系统

奇异值分解实验:图像压缩与推荐系统_奇异值分解_91

您看,我正在看电影,右边会有一个推荐列表。

TA这个是根据什么推荐呢?

可能是我根据的观看记录,这里我们以评分为判断依据吧,简单起见。

奇异值分解实验:图像压缩与推荐系统_工程应用_92


上图一共 奇异值分解实验:图像压缩与推荐系统_奇异值分解_93 位大佬,一共有 奇异值分解实验:图像压缩与推荐系统_工程应用_94 部电影,电影评分是 奇异值分解实验:图像压缩与推荐系统_特征值_95奇异值分解实验:图像压缩与推荐系统_特征值_96

现在【冯八】大佬又来看电影了,我们应该推荐什么给【冯八】呢?

因为每一部电影都有分类的,我们可以在一个类别里面给【冯八】挑:

  • 科幻:变形金刚、钢铁侠、流浪地球
  • 喜剧:喜剧之王、功夫、少林足球
  • 赌博:赌侠、赌神、赌圣

【冯八】给赌博片的打分普遍很高(赌圣5分、赌神4分),所以应该推荐赌博片里没看过的赌侠。

不过计算机理解不了这个影片分类,所以我们可以使用降维,使得数据投影变成 奇异值分解实验:图像压缩与推荐系统_奇异值分解

奇异值分解实验:图像压缩与推荐系统_奇异值分解_98

介绍一下,推荐系统的流程:

  • 前置准备,把所有用户的评分数据放入到一个矩阵里。

奇异值分解实验:图像压缩与推荐系统_特征值_99

  1. 现在为【冯八】推送服务,寻找【冯八】未打分的电影 — 在这个矩阵的第九行里寻找等于 奇异值分解实验:图像压缩与推荐系统_特征值_100

奇异值分解实验:图像压缩与推荐系统_特征值_101

  1. 预测【冯八】会给那些未打分的电影,打多少分。

    科幻类(黄色)打 奇异值分解实验:图像压缩与推荐系统_奇异值分解_102 分,喜剧类(蓝色)打 奇异值分解实验:图像压缩与推荐系统_奇异值分解_102 分,赌博片(红色)打 奇异值分解实验:图像压缩与推荐系统_奇异值分解_104

这里的分类,其实是降维操作 — 使用 奇异值分解实验:图像压缩与推荐系统_工程应用_105 算法将高维投影低维。奇异值分解实验:图像压缩与推荐系统_工程应用_105 算法可参考《特征值分解实验:人脸识别与PageRank网页排序》。不过,这里面的 奇异值分解实验:图像压缩与推荐系统_工程应用_105 算法采用的是特征值分解(EVD),这篇文章是奇异值分解(SVD),所以我们还是用奇异值分解包装的 奇异值分解实验:图像压缩与推荐系统_工程应用_105

奇异向量也可以构造 降维矩阵,因为我们现在是对行降维,协方差矩阵维是 奇异值分解实验:图像压缩与推荐系统_奇异值分解_109,做奇异值分解时,左奇异矩阵 奇异值分解实验:图像压缩与推荐系统_奇异值分解_44 就是矩阵奇异值分解实验:图像压缩与推荐系统_工程应用_111

所以,我们使用 奇异值分解实验:图像压缩与推荐系统_工程应用_112 来构造降维矩阵,那降维矩阵就是从 左奇异矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_44 中截取:奇异值分解实验:图像压缩与推荐系统_工程应用_114

因为 左奇异矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_44 是列向量横向堆叠而成的,所以要转置一下。而后将 数据矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_41 投影到 奇异值分解实验:图像压缩与推荐系统_奇异值分解_44 所代表的低维空间里,得到矩阵奇异值分解实验:图像压缩与推荐系统_工程应用_118

如果是对行降维,协方差矩阵维是 奇异值分解实验:图像压缩与推荐系统_奇异值分解_109,降维矩阵从左奇异矩阵奇异值分解实验:图像压缩与推荐系统_奇异值分解_44 中截取:奇异值分解实验:图像压缩与推荐系统_工程应用_114

如果是对列降维,协方差矩阵维是 奇异值分解实验:图像压缩与推荐系统_工程应用_122,降维矩阵从右奇异矩阵奇异值分解实验:图像压缩与推荐系统_工程应用_52 中截取:奇异值分解实验:图像压缩与推荐系统_特征值_124

在低维空间中,计算出待预测电影与其他电影的相似度。计算相似度,采用相似度算法接口,可参考《向量实验:相似度算法》。

而后,逐一将已评分电影的分数 奇异值分解实验:图像压缩与推荐系统_工程应用_125

  1. 根据预测评分的大小排序,就前 奇异值分解实验:图像压缩与推荐系统_特征值_126

    完整代码:
import numpy as np
 
def load_dataSet():  # “用户-电影”矩阵 ,行表示用户的评分 ,列表示电影
    return np.mat([   [ 5, 4, 5, 4, 0, 0, 0, 0, 1],
                       [0, 0, 0, 0, 5, 4, 1, 4, 0],
                       [0, 0, 0, 0, 0, 0, 0, 5, 4],
                       [3, 3, 5, 0, 5, 0, 0, 0, 0],
                       [0, 0, 0, 0, 1, 1, 0, 0, 0],
                       [1, 2, 3, 0, 0, 0, 0, 0, 0],
                       [2, 0, 0, 5, 4, 5, 0, 0, 0],
                       [0, 0, 5, 0, 0, 0, 0, 1, 0],
                       [1, 0, 2, 0, 1, 0, 0, 4, 5],
                       [0, 5, 0, 0, 0, 0, 0, 1, 0],
                       [4, 4, 4, 0, 0, 0, 0, 1, 2]])
 
def cosSim(inA, inB): 
    return 0.5 + 0.5 * (float(inA.T * inB) / (np.linalg.norm(inA) * np.linalg.norm(inB)))  # 归一化到[0,1],配合分数的归一化
    
def scorePredict(dataMat, xformedItems, user_id, unrated_idx):
    rateTotal = 0.0   # 预测分数
    simTotal = 0.0    # 总相似度(权重)
    n = dataMat.shape[1]        # 获取电影个数
    for i in range(n):         # 遍历所有电影
        userRating = dataMat[user_id, i]  # 针对该用户,拿到一个电影得分   [1, 0, 2, 0, 1, 0, 0, 4, 5],
        if userRating == 0 :   # 跳过未评分项
            continue
        similarity = cosSim(xformedItems[:, unrated_idx], xformedItems[:, i])   # 求余弦相似度
        print( 'the movie_%d and movie_%d similarity is: %f' % (unrated_idx, i, similarity))  
        rateTotal += similarity * userRating  # 预测分数 = 相似度 * 已评分数
        simTotal  += similarity               # 相似度求和
        
    return rateTotal / simTotal  # 评分归一化:使得评分值在0-5之间    
    
    
def  recommed(dataMat,user_id,N=3):
    # 1. 找出该用户未评分电影     
    unratedItems = np.nonzero(dataMat[user_id, :]==0)[1]  # “==0”操作将0置为1,将非0置为0  [0, 1, 0, 1, 0, 1, 1, 0, 0]
    print("-------- The user -------\n",np.around(dataMat[user_id, :], decimals=3)) 
    print("-------- unratedItems -------\n",np.around(unratedItems , decimals=3)) 
    
    # 2.预测评分
    # 2.1. 降维(提取电影主题)
    U, Sigma, VT = np.linalg.svd(dataMat)
    # U*U.T = E ?若为E证明U为正交矩阵,其列向量已经单位正交化,就不用像EVD降维那样,还要自己单位化
    print("----- U*U.T = E ? -----\n",np.around(U*U.T, decimals=0)) 
    
    # 2.2 自动收缩最适合的k
    k = 0  
    for i in range(len(Sigma)):
        if (np.linalg.norm(Sigma[:i + 1]) / np.linalg.norm(Sigma)) > 0.9:  
            k = i + 1
            break   #刚好找到满足条件的k,退出循环
    
    # 2.3 截取U,得到降维矩阵
    red_U = U[:, :k]
    
    # 2.4 降维
    xformedItems = red_U.T * dataMat
    print("xformedItems shape:",xformedItems.shape)  # (3, 9)
    print("----- xformedItems -----\n",np.around(xformedItems, decimals=2)) 
    
    # 2.5 对未评分电影逐一进行分数预测
    movScores = [] # 存储预测到的分数
    for unrated_idx in unratedItems:  # 遍历所有未评分项的索引,逐项预测
        print ("-------- predict movie_%d -------" % (unrated_idx))
        score = scorePredict(dataMat, xformedItems, user_id, unrated_idx)     # 预测当前未评分项的分数
        movScores.append((unrated_idx, score))  # 以元组方式堆叠到movScores
    print("-------- movScores -------\n",np.around(movScores, decimals=3))
    
    # 3.按照预测分数从大到小排序,并返回前N大分数对应的电影
    return sorted(movScores, key=lambda tmp: tmp[1], reverse=True)[:N]  
    
    
if __name__ == "__main__":
    # 1.加载数据集
    dataMat = load_dataSet()
    print("dataMat shape:",dataMat.shape)
    print("-------- dataMat -------\n",np.around(dataMat, decimals=3))
    
    # 2.输入一个用户编号,给他推荐N部电影
    user_id = 8
    N = 4
    recommed_items = recommed(dataMat,user_id,N)
    print("---the recommendation of our system for user_%d are as follows---"%(user_id))
    print(np.around(recommed_items, decimals=3))
    print("done!!!!")

就我们设计的推荐系统,可能面临的一些问题:

  • 可能会有上亿用户,数据矩阵规模很大,矩阵分解会很耗时间;
    解决:因为这个矩阵在一段时间之内变换不大,所以一般一天计算一次就好。
  • 电影有成千上万部,需要多次计算相似度,也很耗时间;
    解决:提前计算各个电影直接的相似度,需要的时候调用即可,不用计算。
  • 很多用户都没有给电影打分的习惯,所以矩阵爆奇异值分解实验:图像压缩与推荐系统_特征值_127,会影响推荐效果。
    解决:胡歌,请您来打分~