清华大学研究生高等数值分析计算实验奇异值分解SVD以及图像压缩matlab源程序代码

第1部分方法介绍奇异值分解（SVD）定理：设，则存在正交矩阵和，使得其中，而且，称为的奇异值，的第列称为的左奇异向量，的第列称为的右奇异向量。注：不失一般性，可以假设，（对于的情况，可以先对转置，然后进行分解，最后对所得的分解式进行转置，就可以得到原来的分解式）方法1：传统的算法主要思想：设，先将二对角化，即构造正交矩阵和使得其中然后，对三角矩阵进行带位移的对称迭代得到。当某个时，具有形状，此时可以将的奇异值问题分解为两个低阶二对角阵的奇异值分解问题；而当某个时，可以适当选取变换，使得第行元素全为零的二对角阵，因此，此时也可以将约化为两个低阶二对角阵的奇异值分解问题。在实际计算时，当或者（这里是一个略大于机器精度的正数）时，就将或者视作零，就可以将分解为两个低阶二对角阵的奇异值分解问题。主要步骤：（1）输入及允许误差（2）计算变换，，使得其中，（3）收敛性检验：（i）将所有满足的置零；（ii）如果，则输出有关信息结束；否则，，确定正整数，使得，，；（iii）如果存在满足使得，则，，，，，转步（iv）；否则转步（4）（iv）确定和使这也相应于所以可以直接调用变换算法得到，这相当于（v）如果，则，，，转步（iv），否则转步（i）（4）构造正交阵和，使仍为上双对角阵，其中，得，，然后转步（3）。方法2：分而治之方法先将矩阵二对角化得到二对角阵，然后将求的问题分为两个子问题。先将矩阵写成如下形式：其中，为上二对角矩阵，是相应维数的向量中的第个单位向量，并且一般我们取。现在假如我们已经求得了，的如下：，并且令为的最后一行，为的第一行。那么将这些带入，可以得到：可以看出中间的矩阵形式上很简单，除了对角元素和第行上的元素，其它元素都为；这里先把它的的计算问题放到最后，而假定已知其为：，将它带入上式，就可以得到的分解式：其中；而计算，的时也可递归地应用这种办法，至到子问题足够得小。现在来讨论中间矩阵的问题，注意到通过排序，将第行和第列移到第1行和第1列，得到矩阵：其中是矩阵，的对角元，而是中间矩阵的第行元素，其中为原位置上的元素；我们继续对矩阵进行排序，并且为了方便定义：，使得排序后，，进而可以方便的求解原矩阵的分解。方法3：方法主要思想：传统的算法和分而治之方法都是基于先将矩阵二对角化的方法，而方法是完全不同的方法，它通过一系列平面旋转（也称为变换）；最终使得矩阵收敛于一个对角矩阵。其中每一个设计成，使得对于一对选定的指标，有如下式子成立：最终收敛时，令，有，且具有互相正交的列向量，可以将写成，其中为对角矩阵，为正交矩阵；从而得到原矩阵的分解式。主要步骤：令是最初的矩阵，是经过一次旋转过得到的矩阵，其中为对角矩阵，是列向量范数全为1的矩阵。可以知道对进行单边的变换相当于对矩阵进行双边的变换。最终得到的矩阵的列向量的范数就是所要求的奇异值，并且最终归一化以后的矩阵的列向量就是左奇异向量。第2部分实验计算针对以上三种算法用进行数值模拟，三个算法源自工具箱，分别为（传统的SVD算法），（分而治之）和（方法）。实验中找10幅彩色图片（尺寸为），并将其灰度化，转化为数据矩阵，计算3种算法的计算时间和最小奇异值误差，以及左右奇异向量误差，并得到了相应的对比曲线图。左右奇异向量的误差通过和的矩阵2-范数的来评估。最小特征值误差，使用默认自带的SVD算法得到的结果作为参考标准（实际中大多数情况也是无法知道确切的值）。1.计算时间的比较：比较结果可以看出，方法的运算时间最长，传统的算法和分而治之方法要快得多，其中分而治之方法是最快的，所需时间不到方法的。2.最小奇异值误差比较结果可以看出，3种算法对最小奇异值误差的影响差别不是很大，3种算法的误差都在以内。3.左奇异向量误差比较结果可以看出，方法的左奇异向量误差最大，传统的算法和分而治之方法要小得多，其中分而治之方法是最小的。4.右奇异向量误差比较结果可以看出，方法和传统的算法的右奇异向量误差较大，而分而治之方法要小得多。小结：传统的SVD算法，比较成熟，但是性能和精度上都不是很理想；方法的效率不高，计算量大，所需的计算时间长；分而治之方法是一种快速算法，但相对精度不一定可靠，主要是由于其中的好多细节都因为技术没公开而缺少文献资料。实际应用中，可以根据需要选择特定的SVD算法。（注：要成功运行以上3种算法，Matlab程序中需要装有NAG工具箱）第3部分图像压缩找5幅彩色图片（尺寸为），并将其转化为数据矩阵，对这些矩阵做分解，并利用公式：这里是奇异值，，为对应的奇异向量，，且。选择，构造，并将其转化为图像输出。本次试验取给出压缩后的图像，并和原图进行比较。彩色图片1：原始彩色图片时时时彩色图片2：原始彩色图片时时时彩色图片3：原始彩色图片时时时彩色图片4：原始彩色图片时时时彩色图片5：原始彩色图片时时时对比分析，可以看出时，图像很不清晰，只能看个大概；时，图像已经比较清晰了，但仍有些模糊；时，图像已经和原图差不多了，肉眼已经很难看出与原图的差别了。程序文件说明svdk.m图像压缩中进行svd分解main.m图像压缩的主程序main2.m比较3种算法的主程序注：运行main2.m程序，Matlab程序中需要装有NAG工具箱。svdk.mfunctionB=svdk(A,k)A=double(A);[U1,S1,V1]=svds(A(:,:,1),k);a1=U1*S1*V1';[U2,S2,V2]=svds(A(:,:,2),k);a2=U2*S2*V2';[U3,S3,V3]=svds(A(:,:,3),k);a3=U3*S3*V3';B=cat(3,a1,a2,a3);B=uint8(B);endmain.mclearall;closeall;clc;N=5;str='.jpg';fori=1:Nch=num2str(i);filename=strcat(ch,str);A=imread(filename);k=20;B=svdk(A,k);ch1=num2str(k);filename1=strcat(ch,'-',ch1,str);imwrite(B,filename1);k=80;B=svdk(A,k);ch1=num2str(k);filename1=strcat(ch,'-',ch1,str);imwrite(B,filename1);k=300;B=svdk(A,k);ch1=num2str(k);filename1=strcat(ch,'-',ch1,str);imwrite(B,filename1);endmain2.mclearall;closeall;clc;N=10;str='.jpg';A0=zeros(N,1);A1=zeros(N,3);A2=zeros(N,3);A3=zeros(N,3);A4=zeros(N,3);fori=1:Nch=num2str(i);filename=strcat(ch,str);A=imread(filename);A=rgb2gray(A);A=A';A=double(A);[m,n]=size(A);I1=ones(m,1);I1=diag(I1);I2=ones(n,1);I2=diag(I2);R=rank(A);A0(i)=R;[U,S,V]=svd(A);S=diag(S);t=cputime;[~,S1,U1,V1,~,~]=f08kb('A','A',A);t=cputime-t;A4(i,1)=t;t=cputime;[~,S2,U2,V2,~]=f08kd('A',A);t=cputime-t;A4(i,2)=t;t=cputime;[U3,S3,V3,~,~]=f08kj('G','U','V',A,int32(0),zeros(n),zeros(m+n,1));t=cputime-t;A4(i,3)=t;A1(i,1)=abs(S(R)-S1(R));A1(i,2)=abs(S(R)-S2(R));A1(i,3)=abs(S(R)-S3(R));A2(i,1)=norm(U1'*U1-I1);A2(i,2)=norm(U2'*U2-I1);A2(i,3)=norm(U3'*U3-I2);A3(i,1)=norm(V1'*V1-I2);A3(i,2)=norm(V2'*V2-I2);A3(i,3)=norm(V3'*V3-I2);endTT=1:N;plot(TT,A1(:,1),TT,A1(:,2),TT,A1(:,3))legend('传统的SVD算法','分而治之','Jacobi');title('最小奇异值误差');plot(TT,A2(:,1),TT,A2(:,2),TT,A2(:,3))legend('传统的SVD算法','分而治之