一、引入

主成分分析PCA、潜在语义分析都会用到SVD
不要求A矩阵是方阵，SVD是线性代数中相似对角化的延伸
任意mn的矩阵都可以用三个矩阵相乘的形式表示
分别是m阶正交矩阵、由降序排列的非负的对角线元素构成的mn矩形对角矩阵、n阶正交矩阵
矩阵奇异值分解一定存在但不唯一
SVD可以看作矩阵压缩数据的一种方法，这种近似是在平方损失意义下的最优近似

二、SVD的定义、性质

定义

A=UΣV^T
满足以下条件
UU^T=E;
VV^T=E;
Σ=diag(σ₁,σ₂,…σ_p),其中σ₁>=σ₂>=…>=σ_p>=0；p=min{m,n}

例题

奇异值分解一定存在

设定m>=n,分三步完成证明
1、确定v和Σ
①Σ
A^TA的特征值都是非负的||Ax||²=x^TA^TAx=λx^Tx=λ||x||²
所以λ=|Ax||²/||x||²>=0
假设正交矩阵V的列的排列使对应的特征值降序排列：λ₁>=λ₂>=λ₃>=…>=λ_n
计算矩阵A的奇异值σ_j=√λj, j=12,…n
说明：A^TA的特征值与A的特征值是平方关系
A^TA=(UΣV^T)^T(UΣV^T)=V(Σ^TΣ)V^T;
AA^T=(UΣV^T)(UΣV^T)^T=U(ΣΣ^T)U^T;
V的列向量是A^TA的特征向量，U的列向量是AA^T的特征向量，Σ的奇异值是A^TA和AA^T的特征值的平方根，A^TA和AA^T的特征值相同

R(A)=r,R(A^TA)=r,由于A^TA是对称矩阵，它的秩等于正的特征值的个数
所以剩下的n-r个特征值为0。所以σ同理

Σ= ②v
V₁=[v₁,v₂,…v_r];V₂=[v_r+1,…v_n]
其中V₁对应的是A^TA的正特征值对用的特征向量， 其中V₂对应的是A^TA的0特征值对用的特征向量V=[V₁,V₂]
V₂
A^TAx=0(V₂的列向量构成了A^TA的零空间N(A^TA)=N(A),所以V₂的列向量构成A的零空间的一组标准正交基)
Ax=0（正交矩阵的转置乘以正交矩阵等于单位矩阵）
③U
u_j=1/σ_jAv_j,j=1,2…r
U₁=[u₁,u₂,…u_r]
则AV₁=U₁Σ₁
U₂与V~2同理

紧奇异值分解和截断奇异值分解

紧奇异值分解（无损压缩，与原始矩阵秩相同）
截断奇异值分解（有损压缩，小于原始矩阵的秩）：满足了秩的要求以后其余元素都变成零

几何解释

被分成的三个矩阵可以解释为，一个坐标轴的旋转或反射变换、一个坐标轴的缩放变换、另一个坐标轴的旋转或反射变换。

三、SVD算法

计算过程

①计算A^TA的特征值和特征向量，特征值开方从大到小排序即为Σ
注意：因为A不是方阵，所以在构造的时候，要将格式修改为m*n的形式，缺的位置补零。
②求n阶正交矩阵V:特征向量单位化
③求m阶正交矩阵U:
u_j=1/σ_jAv_j (j=1,2,3…,r) U₁=[u₁,u₂,u₃,…u_r]
求A^T的零空间的一组标准正交基{u_r+1,u_r+2,…u_m}
（A^Tx=0,求出特征向量以后记得标准化）
U=[U₁,U₂]
④得到奇异值分解

四、SVD与矩阵近似

弗罗贝尼乌斯范数：是向量L₂范数的直接推广，对应机器学习中的平方损失函数

矩阵的外积展开式

若A的秩为n，A_k的秩为k, 且A_k是秩为k的矩阵中在弗罗贝尼乌斯范数意义下A的最优近似矩阵。那么A_k就是A的截断奇异值分解。
通常奇异值σ_i递减的很快，所以k取很小值时，A_K也可以对A有很好的近似。

五、python实现

import numpy as np
a = np.random.randint(-10,10,(4, 3)).astype(float)
print(a)
print("-----------------")
u, sigma, vT = np.linalg.svd(a)
print(u)
print("-----------------")
print(sigma)
print("-----------------")
print(vT)
print("-----------------")# 将sigma 转成矩阵
SigmaMat = np.zeros((4,3))
SigmaMat[:3, :3] = np.diag(sigma)
print(SigmaMat)
print("------验证-------")
a_ = np.dot(u, np.dot(SigmaMat, vT))
print(a_)

六、应用

推荐算法

 import numpy as npimport randomclass SVD:def __init__(self,mat,K=20):self.mat=np.array(mat)self.K=Kself.bi={}self.bu={}self.qi={}self.pu={}self.avg=np.mean(self.mat[:,2])for i in range(self.mat.shape[0]):uid=self.mat[i,0]iid=self.mat[i,1]self.bi.setdefault(iid,0)self.bu.setdefault(uid,0)self.qi.setdefault(iid,np.random.random((self.K,1))/10*np.sqrt(self.K))self.pu.setdefault(uid,np.random.random((self.K,1))/10*np.sqrt(self.K))
def predict(self,uid,iid):  #预测评分的函数#setdefault的作用是当该用户或者物品未出现过时，新建它的bi,bu,qi,pu，并设置初始值为0self.bi.setdefault(iid,0)self.bu.setdefault(uid,0)self.qi.setdefault(iid,np.zeros((self.K,1)))self.pu.setdefault(uid,np.zeros((self.K,1)))rating=self.avg+self.bi[iid]+self.bu[uid]+np.sum(self.qi[iid]*self.pu[uid]) #预测评分公式#由于评分范围在1到5，所以当分数大于5或小于1时，返回5,1.if rating>5:rating=5if rating<1:rating=1return ratingdef train(self,steps=30,gamma=0.04,Lambda=0.15):    #训练函数，step为迭代次数。print('train data size',self.mat.shape)for step in range(steps):print('step',step+1,'is running')KK=np.random.permutation(self.mat.shape[0]) #随机梯度下降算法，kk为对矩阵进行随机洗牌rmse=0.0for i in range(self.mat.shape[0]):j=KK[i]uid=self.mat[j,0]iid=self.mat[j,1]rating=self.mat[j,2]eui=rating-self.predict(uid, iid)rmse+=eui**2self.bu[uid]+=gamma*(eui-Lambda*self.bu[uid])  self.bi[iid]+=gamma*(eui-Lambda*self.bi[iid])tmp=self.qi[iid]self.qi[iid]+=gamma*(eui*self.pu[uid]-Lambda*self.qi[iid])self.pu[uid]+=gamma*(eui*tmp-Lambda*self.pu[uid])gamma=0.93*gammaprint('rmse is',np.sqrt(rmse/self.mat.shape[0]))def test(self,test_data):  #gamma以0.93的学习率递减test_data=np.array(test_data)print('test data size',test_data.shape)rmse=0.0for i in range(test_data.shape[0]):uid=test_data[i,0]iid=test_data[i,1]rating=test_data[i,2]eui=rating-self.predict(uid, iid)rmse+=eui**2print('rmse of test data is',np.sqrt(rmse/test_data.shape[0]))def getData():   #获取训练集和测试集的函数import ref=open('C:/Users/xuwei/Desktop/data.txt','r')lines=f.readlines()f.close()data=[]
for line in lines:list=re.split('\t|\n',line)if int(list[2]) !=0:    #提出评分0的数据，这部分是用户评论了但是没有评分的data.append([int(i) for i in list[:3]])random.shuffle(data)train_data=data[:int(len(data)*7/10)]test_data=data[int(len(data)*7/10):]print('load data finished')print('total data ',len(data))
return train_data,test_datatrain_data,test_data=getData()
a=SVD(train_data,30)  
a.train()
a.test(test_data)

代码链接：https://blog.csdn.net/akiyamamio11/article/details/79042688