奇异值分解（三）

最新推荐文章于 2024-05-31 00:21:41 发布

我闻如是

最新推荐文章于 2024-05-31 00:21:41 发布

阅读量2.1k

点赞数

分类专栏：压缩感知文章标签：机器学习矩阵线性代数

本文链接：https://blog.csdn.net/xfsong2012/article/details/121318803

版权

压缩感知专栏收录该内容

6 篇文章 0 订阅

订阅专栏

对下面的文章程序进行了实现。
浅谈张量分解（三）：如何对稀疏矩阵进行奇异值分解？
SVD的应用1：矩阵缺失值补全
有一个5*4的稀疏矩阵A
其中“？”代表缺失值，矩阵A的真实值为在这里插入图片描述
这是交通数据，其中行代表的是5个不同的传感器，列代表4个不同的观测时间间隔，即时间窗，比如15min；矩阵里面的数据代表的是交通流量（辆/15min）。
稀疏矩阵的奇异值分解将分为两部分：缺失元素的初始化填补和迭代的奇异值分解。
1、缺失元素的初始化填补
（1）先需要对缺失的元素进行初始化，最简单的可以全部选择0。
程序:

A_sparse=[0,90,449,517;0,0,412,0;192,0,697,687;185,0,699,657;164,58,0,0];
A_real=[208,90,449,517;104,43,412,411;192,77,697,687;185,115,699,657;...
    164,58,696,599];
pos1=find(A_sparse==0);%找出等于0的项
RMSE=sqrt(sum((A_real(pos1)-A_sparse(pos1)).^2)/length(pos1))

计算出RMSE=368.6497，误差太大。这是因为各行列损失程度不同导致的结果，所以不能直接用0替代缺失值。
（2）可以考虑使用均值填补，求出原矩阵的均值=401.用401替代原来的缺失值，并计算RMSE。
程序：

A_sparse=[0,90,449,517;0,0,412,0;192,0,697,687;185,0,699,657;164,58,0,0];
A_real=[208,90,449,517;104,43,412,411;192,77,697,687;185,115,699,657;...
    164,58,696,599];
pos1=find(A_sparse==0);%找出等于0的项
pos2=find(A_sparse~=0);%找出非0的项
[m,n]=size(A_real);%算出A实际值的行列
mu=401;
error=zeros(m,n);
 for i=1:m
        for j=1:n
            if A_sparse(i,j)==0
                error(i,j)=mu;
            end
        end
 end
 A_hat=error+A_sparse
RMSE=sqrt(sum((A_real(pos1)-A_hat(pos1)).^2)/length(pos1))

输出RMSE=266.1586，误差还是比较大。
（3）为了求取更精确的值，可以用如下公式对估计值进行修正：
在这里插入图片描述
程序：

clear
clc
A_sparse=[0,90,449,517;0,0,412,0;192,0,697,687;185,0,699,657;164,58,0,0];
A_real=[208,90,449,517;104,43,412,411;192,77,697,687;185,115,699,657;...
    164,58,696,599];
pos1=find(A_sparse==0);%矩阵中值为0的索引
pos2=find(A_sparse~=0);%矩阵中值不为0的索引
[m,n]=size(A_real);%求出行和列
mu=mean(A_sparse(pos2));%求出缺失矩阵的均值
b1=0.1*rand(1,m);%行上的修正值
b2=0.1*rand(1,n);%列上的修正值
error=zeros(m,n);%误差
alpha1=0.01;%学习率
for iter=1:1000 %迭代次数
    for i=1:m
        for j=1:n
            A_hat(i,j)=round(mu+b1(1,i)+b2(1,j));%求估计值并取整
            if A_sparse(i,j)~=0
                error(i,j)=A_sparse(i,j)-A_hat(i,j);
            end
        end
    end
    b1=b1+alpha1*sum(error');b2=b2+alpha1*sum(error);
    L1(1,iter)=0.5*sum(sum(error.^2));
end
A_hat
RMSE=sqrt(sum((A_real(pos1)-A_hat(pos1)).^2)/length(pos1))

输出：
A_hat =

48     41   505  510
-44   -51   412  418
219   212   676  681
207   200   664  670
115   107   571  577

RMSE = 110.6504。
虽然RMSE值变小了，但是估计的A值里面有负值。这种结果针对于车流量来说不合理，预估结果不理想。
（4）采用下式进行修补
在这里插入图片描述
程序：

A_sparse=[0,90,449,517;0,0,412,0;192,0,697,687;185,0,699,657;164,58,0,0];
A_real=[208,90,449,517;104,43,412,411;192,77,697,687;185,115,699,657;...
    164,58,696,599];
pos=find(A_sparse==0);
[m,n]=size(A_real);
r=3;
U=0.1*rand(m,r);
V=0.1*rand(n,r);
error=zeros(m,n);
alpha2=0.0001;
for iter=1:800
    A_hat=round(U*V');
    error=A_sparse-A_hat;
    error(pos)=0;
    U_plus=U+alpha2*error*V;
    V_plus=V+alpha2*error'*U;
    U=U_plus;
    V=V_plus;
    L2(1,iter)=0.5*sum(sum(error.^2));
end
A_hat
RMSE=sqrt(sum((A_real(pos)-A_hat(pos)).^2)/length(pos))

输出：
A_hat =

135 72 476 492
116 57 412 410
194 95 693 690
191 91 682 672
154 74 551 545

RMSE =

61.7566
经过修正后，误差已经小了很多。
2、迭代的奇异值分解
在这里插入图片描述
程序：

clear
clc
A_sparse=[0,90,449,517;0,0,412,0;192,0,697,687;185,0,699,657;164,58,0,0];
A_real=[208,90,449,517;104,43,412,411;192,77,697,687;185,115,699,657;...
    164,58,696,599];
A_hat=[401,90,449,517;401,401,412,401;192,401,697,687;185,401,699,657;...
    164,58,401,401];
pos=find(A_sparse==0);
[m,n]=size(A_real);
k=2;
A_f=A_hat;
RMSE=zeros(1);
for iter=1:8000
    [Q,Sigma,P]=svds(A_f,k);
    A_hat=Q*Sigma*P';
    A_f(pos)=A_hat(pos);
    convergence(1,iter)=sum(A_hat(pos).^2);
end
RMSE=sqrt(sum((A_real(pos)-A_f(pos)).^2)/length(pos))

输出RMSE=161.848
用（3）中得到的A估计值。

A_sparse=[0,90,449,517;0,0,412,0;192,0,697,687;185,0,699,657;164,58,0,0];
A_real=[208,90,449,517;104,43,412,411;192,77,697,687;185,115,699,657;...
    164,58,696,599];
A_hat=[48,41,505,510;-44,-51,412,418;219,212,676,681;207,200,664,670;...
    115,107,571,577];
pos=find(A_sparse==0);[m,n]=size(A_real);k=2;
A_f=A_hat;
for iter=1:100
    [Q,Sigma,P]=svds(A_f,k);
    A_hat=Q*Sigma*P';
    A_f(pos)=A_hat(pos);
    convergence(1,iter)=sum(A_hat(pos).^2);
end
RMSE=sqrt(sum((A_real(pos)-A_f(pos)).^2)/length(pos))

RMSE=110.15
用（4）中的A评估值

A_sparse=[0,90,449,517;0,0,412,0;192,0,697,687;185,0,699,657;164,58,0,0];
A_real=[208,90,449,517;104,43,412,411;192,77,697,687;185,115,699,657;...
    164,58,696,599];
A_hat=[139,89,449,517;114,46,412,408;192,75,697,687;186,59,699,657;...
    163,59,599,579];
pos=find(A_sparse==0);[m,n]=size(A_real);k=2;
A_f=A_hat;
for iter=1:500
    [Q,Sigma,P]=svds(A_f,k);
    A_hat=Q*Sigma*P';
    A_f(pos)=A_hat(pos);
    convergence(1,iter)=sum(A_hat(pos).^2);
end
RMSE=sqrt(sum((A_real(pos)-A_f(pos)).^2)/length(pos))