清华大学高等数值计算(李津)实践题目二(SVD计算及图像压缩)(包含matlab代码)

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