奇异值分解实验
- 奇异值分解
- 低秩近似
- 工程应用:图像压缩
- 工程应用:推荐系统
奇异值分解
只有方阵(行数等于列数)才能做特征值分解,非方阵可不可以分解为
这种方式是【奇异值分解】,这种方法大学里并不学。
因为本科的线性代数主要研究方阵(除了线性系统),所以大学里并没有介绍非方阵的奇异值分解(),奇异值分解在数据降维、语义分析、图像等领域都有十分广泛的应用,比如
举个荔枝,演示一下非方阵的分解步骤。
- 求出矩阵 和 :
发现 和
- 因此,分别求出 和
有俩组,
有三组,
- 将 的特征向量横向拼成矩阵 :
- 将 和 的相同特征值开方,拼成矩阵:
- 拼成矩阵:
P.S.
- 将 的特征向量横向拼成一个矩阵:
P.S. 也是正交矩阵,因为TA的列向量俩俩正交,且是单位向量。
- 最后,将 发现结果等于原矩阵,说明 矩阵 可以分解为。
虽然这只是一个 的矩阵的分解过程,但推广到
- 奇异值分解:
, 矩阵 的尺寸是
,矩阵 的尺寸是 ,其列向量称为 矩阵
,矩阵 的尺寸是 ,其对角线上的值称为 矩阵
, 矩阵 的尺寸是 ,其列向量称为 矩阵
低秩近似
我们可以把矩阵看成一种变换,把矩阵乘法当成线性变换时,找出变换矩阵的特征值和特征向量,实际上就是找出变换矩阵的主要变换方向。
也可以说是,特征值和特征向量代表了一个方阵的【固有信息】。
特征值分解是奇异值分解的特例,特征值分解只能分解方阵,奇异值分解可以分解任意形状的矩阵。
因此,奇异值及奇异向量可以说是代表了一个
奇异值越大,代表的信息就越多。
另外,如果我们在奇异值矩阵
很多情况下,前 个奇异值的和就占了全部奇异值之和的 %。
也就是说,前
所以,可以用 最大的
这种理论,就被称为【低秩近似】。
来看一个
奇异值分解:,看中间的矩阵: 。
发现奇异值有 个,分别是 ,我们只用前面俩个奇异值来算一下,占总奇异值的比例:
- %
这个比例很大了,所以我们可以认为前面俩个奇异值,及其对应的左右奇异向量足以代表原来的矩阵。
把这些部分截取下来:
将截取下来的矩阵相乘,就低秩近似原来的数据矩阵,对比俩个矩阵可以发现,俩者数值非常接近。
- 原始矩阵: 近似矩阵:
工程应用:图像压缩
基于奇异值分解的【低秩近似】理论在工程中有广泛的应用,比如图像压缩。
图像本来就是一个矩阵,比如下面的图片:
假设这张图片的尺寸是 ,就需要
就这样一张不大的灰度图,都要将近 万字节()存储,要是某个应用是图片为主的,那您可以想象应用会有多大。
图像压缩主要有俩个好处:
- 存储空间会小很多
- 方便网络传输
我们可以用 奇异值分解 来压缩图像,算法就是【低秩近似】理论。
图像压缩有俩部分:
- 压缩图像
对 图像矩阵做奇异值分解,得到
选取前 大的奇异值(),按照【低秩近似】理论,对 做截取,得到
存储或者传输新的 ,新的矩阵不一定比原来的小,一定要选取一个恰当的
- 图像重构
将
完整代码:
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
哈哈,整整压缩了一个量级()。
显示图片:
工程应用:推荐系统
您看,我正在看电影,右边会有一个推荐列表。
TA这个是根据什么推荐呢?
可能是我根据的观看记录,这里我们以评分为判断依据吧,简单起见。
上图一共 位大佬,一共有 部电影,电影评分是 ,
现在【冯八】大佬又来看电影了,我们应该推荐什么给【冯八】呢?
因为每一部电影都有分类的,我们可以在一个类别里面给【冯八】挑:
- 科幻:变形金刚、钢铁侠、流浪地球
- 喜剧:喜剧之王、功夫、少林足球
- 赌博:赌侠、赌神、赌圣
【冯八】给赌博片的打分普遍很高(赌圣5分、赌神4分),所以应该推荐赌博片里没看过的赌侠。
不过计算机理解不了这个影片分类,所以我们可以使用降维,使得数据投影变成
介绍一下,推荐系统的流程:
- 前置准备,把所有用户的评分数据放入到一个矩阵里。
- 现在为【冯八】推送服务,寻找【冯八】未打分的电影 — 在这个矩阵的第九行里寻找等于
- 预测【冯八】会给那些未打分的电影,打多少分。
科幻类(黄色)打 分,喜剧类(蓝色)打 分,赌博片(红色)打
这里的分类,其实是降维操作 — 使用 算法将高维投影低维。 算法可参考《特征值分解实验:人脸识别与PageRank网页排序》。不过,这里面的 算法采用的是特征值分解(EVD),这篇文章是奇异值分解(SVD),所以我们还是用奇异值分解包装的
奇异向量也可以构造 降维矩阵,因为我们现在是对行降维,协方差矩阵维是 ,做奇异值分解时,左奇异矩阵 就是矩阵
所以,我们使用 来构造降维矩阵,那降维矩阵就是从 左奇异矩阵 中截取:。
因为 左奇异矩阵 是列向量横向堆叠而成的,所以要转置一下。而后将 数据矩阵 投影到 所代表的低维空间里,得到矩阵。
如果是对行降维,协方差矩阵维是 ,降维矩阵从左奇异矩阵 中截取:。
如果是对列降维,协方差矩阵维是 ,降维矩阵从右奇异矩阵 中截取:。
在低维空间中,计算出待预测电影与其他电影的相似度。计算相似度,采用相似度算法接口,可参考《向量实验:相似度算法》。
而后,逐一将已评分电影的分数
- 根据预测评分的大小排序,就前
完整代码:
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!!!!")
就我们设计的推荐系统,可能面临的一些问题:
- 可能会有上亿用户,数据矩阵规模很大,矩阵分解会很耗时间;
解决:因为这个矩阵在一段时间之内变换不大,所以一般一天计算一次就好。 - 电影有成千上万部,需要多次计算相似度,也很耗时间;
解决:提前计算各个电影直接的相似度,需要的时候调用即可,不用计算。 - 很多用户都没有给电影打分的习惯,所以矩阵爆,会影响推荐效果。
解决:胡歌,请您来打分~