数值分析:矩阵奇异值分解及其应用(Numpy实现)
算法的完整实现代码我已经上传到了GitHub仓库:NumericalAnalysis-Python(包括其它数值分析算法),感兴趣的童鞋可以前往查看。
1 奇异值分解(SVD)
1.1 奇异值分解
已知矩阵A∈Rm×nA∈Rm×n, 其奇异值分解为:
其中U∈Rm×mU∈Rm×m,V∈Rn×nV∈Rn×n是正交矩阵,S∈Rm×nS∈Rm×n是对角线矩阵。SS的对角线元素s1,s2,...,smin(m,n)s1,s2,...,smin(m,n)是矩阵的奇异值。
1.2 奇异值分解的求解
而求矩阵的奇异值的算法非常简单,对于实数域下的矩阵AA,我们只需要求ATAATA的特征值和特征向量。其特征向量归一化后即右奇异向量v1,v2,...,vnv1,v2,...,vn,其特征值开根号即对应的奇异值s1,s2,...,smin(m,n)s1,s2,...,smin(m,n)。 然后由等式
依次计算出相应的ui向量的值。
至于特征值的计算,采用 QR 算法,此处不予介绍,这里可以直接调用 np.linalg.eig()
函数实现。以下给出奇异值计算代码实例(此处仅为知识演示,具体的工业级别的奇异值计算算法要复杂得多,参考 Golub 与 Van Loan《矩阵计算》)
import numpy as np
def svd(A):
eigen_values, eigen_vectors = np.linalg.eig(A.T.dot(A))
singular_values = np.sqrt(eigen_values)
#这里奇异值要从大到小排序,特征向量也要随之从大到小排
val_vec = [] #存储奇异值-特征向量对
for i in range(len(eigen_values)):
val_vec.append((singular_values[i], eigen_vectors[:, i]))
val_vec.sort(key = lambda x:-x[0])
singular_values = [ pair[0] for pair in val_vec]
eigen_vectors = [ pair[1] for pair in val_vec]
# 在计算左奇异向量之前,先要对右奇异向量也就是特征向量组成的基正交化
# 不过linalg.eig返回的是已经正交化的,这一步可省略
# 由等式Avi = siui(vi是右奇异向量, ui是左奇异向量)
# 依次计算左奇异向量
U = np.zeros((A.shape[0], A.shape[1]))
for i in range(A.shape[1]):
u = A.dot(eigen_vectors[i])/singular_values[i]
U[:, i] = u
# 给U加上标准正交基去构造R3的基
for i in range(A.shape[1], A.shape[0]):
basis = np.zeros((A.shape[0], 1))
basis[i] = 1
U = np.concatenate([U, basis], axis=1)
eigen_vectors = [vec.reshape(-1, 1) for vec in eigen_vectors]
eigen_vectors = np.concatenate(eigen_vectors, axis=1)
return U, singular_values, eigen_vectors
if __name__ == '__main__':
# 例一:普通矩阵
A = np.array(
[
[0, 1],
[0, -1]
]
)
# 例二:对称矩阵
# A = np.array(
# [
# [0, 1],
# [1, 3/2]
# ]
# )
U, S, V = svd(A)
print("我们实现的算法结果:")
print(U, "\n", S, "\n", V)
print("\n")
print("调用库函数的计算结果:")
# 调用api核对
U2, S2, V2 = np.linalg.svd(A)
print(U2, "\n", S2, "\n", V2)
对普通矩阵(010−1)运行该算法的结果为:
我们实现的算法结果:
[[ 0.70710678 0.70710678]
[-0.70710678 -0.70710678]]
[1.4142135623730951, 0.0]
[[0. 1.]
[1. 0.]]
调用库函数的计算结果:
[[-0.70710678 -0.70710678]
[ 0.70710678 -0.70710678]]
[1.41421356 0. ]
[[-0. -1.]
[-1. 0.]]
可以看到结果基本符合。此处矩阵(010−1)的奇异值为0.41421和0。此处我们发现普通矩阵的奇异值可以为0。
对对称矩阵(01132)运行该算法的结果为:
我们实现的算法结果:
[[-0.4472136 0.89442719]
[-0.89442719 -0.4472136 ]]
[2.0, 0.5]
[[-0.4472136 -0.89442719]
[-0.89442719 0.4472136 ]]
调用库函数的计算结果:
[[-0.4472136 -0.89442719]
[-0.89442719 0.4472136 ]]
[2. 0.5]
[[-0.4472136 -0.89442719]
[ 0.89442719 -0.4472136 ]]
可以看到结果基本符合。此处矩阵(010−1)的奇异值为2和0.5。此处我们发现,对称矩阵的奇异值必为正,不可能为0。
2 奇异值分解的应用
2.1 奇异值的应用1:推荐系统
在推荐系统中,我们常定义用户-评分矩阵,表示用户对商品的打分,这个矩阵我们称为共现矩阵。
而这就迫切地需要我们设计矩阵分解算法,为每一个用户和视频生成一个隐向量,将用户和视频定位到隐向量的表示空间上,并满足距离相近的用户和视频表示兴趣特点接近。
在推荐系统的应用场景下,我们企图使用矩阵分解算法将m×n维的共现矩阵R分解为m×k维的用户矩阵和k∗n维的物品矩阵(的转置)相乘的形式。其中m是用户数量,n是物品数量,k是隐向量。k的大小决定了隐向量表达能力的强弱。k越小,隐向量包含的信息越少,模型的泛化程度越高;反之,k越大,隐向量表达能力越强,泛化程度相应降低。此外,k的取值还与矩阵分解的求解复杂度直接相关。应用中,k的取值要经过试验多次找到一个推荐效果和工程开销的平衡。具体的形式如下图所示:
采用什么方法来进行矩阵分解呢?由矩阵分析的知识可得,特征值分解只能作用于方阵,显然不适合于分解用户-物品矩阵。我们在这里采用矩阵的奇异值分解以得到用户和物品的隐向量。
已知M是矩阵m×n的矩阵,则一定存在一个分解M=Udiag(λ1,λ2,...,λn)VT,其中U是m∗m的正交矩阵,V是n×n的正交矩阵,diag(λ1,λ2,...,λn) 是 m×n的对角阵。 我们取对角阵 diag(λ1,λ2,...,λn)中较大的k个元素做为隐含特征,删除diag(λ1,λ2,...,λn)中的其他维度及U和V中对应的维度。矩阵M被分解为M=Um∗kdiag(λ1,λ2,...,λk)VTk∗n,至此完成了隐向量维度为k的矩阵分解。 如果我们调用np.lialg.svd()
函数接口,那我们可以将奇异值分解表述如下:
import numpy as np
if __name__ == '__main__':
M = np.array(
[
[0, 4.5, 2.0, 0],
[4.0, 0, 3.5, 0],
[0, 5.0, 0, 2.0],
[0, 3.5, 4.0, 1.0]
]
)
U, S, V_T = np.linalg.svd(M)
k = 2 # 取前2个奇异值对应的隐向量
# 分别打印物品向量和用户向量
Vec_user, Vec_item = U[:,:k], V_T[:k, :].T
print(Vec_user, "\n\n", Vec_item)
该算法对运行结果为:
[[-0.55043774 0.1361732 ]
[-0.26216705 -0.86775439]
[-0.52483774 0.4552962 ]
[-0.5939967 -0.1454804 ]]
[[-0.12135946 -0.63908086]
[-0.83093848 0.43821815]
[-0.50855715 -0.61619448]
[-0.19021762 0.14087178]]
可以看到我们由共现矩阵成功得到了用户向量和物品向量。
然而,运在推荐系统中的传统奇异值分解存在两点重大的缺陷:
- 奇异值分解要求原始共现矩阵是稠密的,而互联网场景下用户非常少,用 户-物品的共现矩阵非常系数。如果使用 SVD,就必须对缺失的元素值进行填充。
- 传统奇异值分解的计算复杂度达到了O(mn2)的级别,这对于商品数量动 辄上百万,用户数量往往上千万的互联网场景来说根本不可接受。 所以,传统奇异值分解不适用于解决大规模稀疏矩阵的矩阵分解。因此,梯度下降法成为了矩阵分解的主要方法。这部分内容我们会在推荐系统专栏中进行讲解。
2.2 奇异值的应用2:矩阵的低秩近似和数据降维
将矩阵的奇异值分解形式M=USVT中的对角阵进一步写成多个子矩阵的和,我们有:
注意,这里u1和v1是做外积,运算得到一个矩阵。 也就是说,m×n的矩阵A可以写成秩为1的矩阵和,即:
我们将这个性质称为 SVD 的低秩近似性质。
在介绍 SVD 的底秩近似的应用前,我们先介绍数据降维的思想。降维的思想是将数据投影到低维空间,假设a1,a2...,an都是m维向量(在数据科学的应用中, 一般m远小于n,想想为什么)。降维的目标是使用n个p维的向量替换原本的n个m维的向量,其中新向量的维度p<m,同时最小化该过程引入的误差。
那么 SVD 其实天然可以用于降维。我们定义矩阵A的秩p近似,将矩阵A的奇异值分解保留前p项,即:
也就是其低秩近似形式保留前p项,
这个式子也可以看做A的最优最小二乘近似形式,即:
这里,B的大小和A一样,B∈Rm×n(但是rank(B)⩽p),F指F范数。这里的F范数可以推广到任意的酉不变范数||A−B||U,不过在常规的使用中,大家就使用F范数就够了。
矩阵最优近似是有着几何解释的。空间<u1,...,up>由左奇异向量u1,...,up长成,这是对于a1,...,an的p维子空间在最小二乘意义上的最优近似,A的列ai在该空间上的正交投影对应Ap的列。换句话讲,一组向量 a1,a2,...,an找到其最优的最小二乘p维子空间的投影(最小二乘后面会介绍,这里暂时理解不了也没关系)就是矩阵最优的秩p近似矩阵Ap。
比如,我们要找到最优的一维子空间拟合数据向量(3,2)T,(2,4)T,(−2,−1)T,(−3,−5)T。 4 个向量近似指向相同的一维子空间,我们想找出这个子空间,该空间能够使向 量投影到子空间的平方误差和最小。然后我们找出投影向量,投影向量组成的矩阵就是我们要求的近似矩阵Ap。
如下图所示:
算法如下:
import numpy as np
from sklearn.decomposition import PCA
def approximation(A, p):
U, s, V_T = np.linalg.svd(A)
B = np.zeros(A.shape)
for i in range(p):
B += s[i]*U[:,i].reshape(-1, 1).dot(V_T[i, :].reshape(1, -1))
return B
if __name__ == '__main__':
# 例一:
# A = np.array(
# [
# [0, 1],
# [1, 3/2],
# ]
# )
# 例二:
A = np.array(
[
[3, 2, -2, -3],
[2, 4, -1, -5]
]
)
# p为近似矩阵的秩,秩p<=r
p = 1
B = approximation(A, p)
print(B)
#最终得到的矩阵秩
print(np.linalg.matrix_rank(B))
(注意,numpy 内置的 SVD 函数返回的是VT而不是V,我就在这儿犯过错。。。 导致后面求出来的近似矩阵不对/(ㄒoㄒ)/~~)
最终对例一矩阵(01132)运行算法的结果如下:
[[0.4 0.8]
[0.8 1.6]]
1
该矩阵的四个列向量对应原始数据向量的投影向量。可以看到这四个向量线性相关, 且最终得到的矩阵的秩为1。
最终对例二矩阵(32−2−324−1−5)运行算法的结果如下:
[[ 1.99120445 2.59641446 -1.16885153 -3.41876737]
[ 2.73456571 3.56571418 -1.60521001 -4.69506988]]
1
该矩阵的四个列向量对应原始数据向量的投影向量。可以看到这四个向量线性相关, 且最终得到的矩阵的秩也为1。
也就是说我们的算法对这两个矩阵都达到了我们低秩近似的效果。因为降维后这两个矩阵的四个向量同属于一个一维子空间,我们只需要一个维度就可以区分这四个向量了,因此我们达到了数据降维的效果。
2.3 奇异值的应用3:压缩
矩阵的奇异值分解可以用于压缩矩阵信息。我们注意到矩阵的展开式
中,每一项使用两个向量ui,vi,以及一个数字si定义。如果A是一个n×n矩阵,我们可以尝试矩阵A的有损压缩,及扔掉求和后面的几项,它们具有较小的si,也就是说对数据的存储而言显得“无关紧要”。就这样,我们可以保留前p项,将矩阵的p秩近似做为矩阵的压缩结果。p越多则近似矩阵对矩阵的近似程度越高,压缩程度越低;p越少则近似矩阵对矩阵的近似程度越低,压缩程度越高。每一项包括ui向量、vi向量和一个数字si,总共需要2n+1个数字保存或者传输。例如,当n=8时,矩阵由64个图片定义。但是我们可以传输或者保存矩阵的第一项展开,仅仅使用2n+1=17个数字。如果大量信息可以由第一项捕捉, 例如,当第一个奇异值比其他的奇异值大得多的时候,以这种方式处理可能节省75%的空间。
算法如下:
import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams["font.sans-serif"] = [u"SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
def approximation(A, p):
B = np.zeros(A.shape)
for c in range(A.shape[2]):
U, s, V_T = np.linalg.svd(A[:, :, c])
for i in range(p):
B[:, :, c] += s[i] * \
U[:, i].reshape(-1, 1).dot(V_T[i, :].reshape(1, -1))
return B
if __name__ == '__main__':
img = cv.imread(
"chapter12.特征值和奇异值/12.4.SVD的应用/12.4.3.图像压缩/img.jpeg")
# 将OpenCV采用的BGR格式转换到Matplotlib采用的RGB格式
img = cv.cvtColor(img, cv.COLOR_BGR2RGB)
# 图像必须归一化到[0 - 1]范围
img = img.astype(np.float32) / 255.0
img_output = img.copy()
# p为近似矩阵的秩,秩p<=r,p越大图像压缩程度越小,越清晰
p = 50
img_output = approximation(img, p)
fig, axs = plt.subplots(1, 2)
axs[0].imshow(np.clip(img, 0, 1))
axs[0].set_title(u"原图")
axs[1].imshow(np.clip(img_output, 0, 1))
axs[1].set_title(u"压缩后的图")
plt.savefig(
"chapter12.特征值和奇异值/12.4.SVD的应用/12.4.3.图像压缩/result.png",
bbox_inches="tight")
plt.show()
最终图片的压缩效果:
3 知名程序库和源码阅读建议
SVD 算法有很多优秀的开源甚至分布式的实现,这里推荐几个项目:
3.1 Gensim
Gensim 是一个采用 Python 和 Cpython 实现的自然语言库,提供了很多统计自然语言处理算法的实现,也包括我们这里提到的 SVD 算法。
文档地址:https://radimrehurek.com/gensim/
源码地址:https://github.com/RaRe-Technologies/gensim.git
3.2 Spark-MLlib
Spark 除了包含 GraphX,它还包括了机器学习库 MLlib,其中就有奇异值分解的分布式实现。
文档地址:https://spark.apache.org/mllib/
源码地址:https://github.com/apache/spark
参考
- [1] Timothy sauer. 数值分析(第2版)[M].机械工业出版社, 2018.
- [2] Golub, Van Loan. 矩阵计算[M]. 人民邮电出版社, 2020.
- [3] 王喆. 深度学习推荐系统[M]. 电子工业出版社, 2020.
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~