清华大学 李津老师 数值分析第二次实验作业

第二次计算实验 SVD 及其应用 一、方法 从已有函数库或软件包(比如 Matlab 库函数及其工具箱，Lapack，Linpack，Propack 等) 选择至少两个用不同方法计算 SVD 的程序，解释算法的思想和主要步骤。 这里，我们主要分析 Matlab 中的 svd()函数和 Propack 中的 lansvd()函数。(lansvd()函数 程序见附件 lansvd.m，svd()由于是 Matlab 中封装好函数无法提取，所以这里不在附附件)。 通过阅读相关文献及程序， 我们发现这两种算法都是通过将原矩阵上双对角化， 再对此 上双对角阵进行奇异值分解，他们主要的区别在于如何将原矩阵上双对角化。 我们令带处理矩阵设为 A，则奇异值分解可描述为下式 0 00 T AVU 1.Matlab 中的 svd()算法描述： (1)先将 m n A (mn)上双对角化，这里使用 Householder 方法，其具体过程如下所示 将 m n A 分块为 11nm AvA 其中 1 v为 m 阶列向量，先计算 m 阶 Householder 变换使得 1 11 111 (,) m PveR eR 并且形成： 11 11 11 1 01 T n mn u APA Am 再计算 n1 阶 Householder 变换 1 H使得 1 112 121 (,) n n H ueR eR 并形成： 12 1 11 1 ,0 10 0 0 T mn A A H 然后对 11mn A 再循环以上运算，共需循环 n-1 次进行后，即可得到上双对角阵 B，以及 U,V 12 23 000 000 000 000 0000 n n B 12 1 .nUP PP，121 1 .nVH HH 1 0 0 i i i I P P ， 1 0 0 i iT i I V H (2)对 B 进行奇异值分解，我们令 T TB B 对 T 进行带 Wilkinson 位移的对称 QR 迭代即可得到 B 的奇异值分解， 进而得到 A 的奇异值 分解。 2.Propack 中 lansvd()算法 (1)使用 lanczos 方法将 m n A 上双对角化 1. 初始化向量 0 m p ， 10 p， 101 up， 0 0v 2. for i=1,2,n 1 111 , , T iiii iiiii iiii iiiii rA uv rvr pAvu pup end 这样迭代 k 步以后，A有如下分解形式： 1nnn AVUB 其中 k V和 1k U 为算法 B.1 中的 Lanczos 向量组合而成，每一列都正交，且已归一化。而 n B为如下形式的二对角矩阵： 1 22 3 1 n n n B 再对 n B进行奇异值分解，其方法与 svd()中对上双对角阵的奇异值分解方法相同，这里 就不再赘述了。 二、实际计算 生成十个不同的(最好属不同类型或有不同性质的)的m n 矩阵，这里,100m n ， 用你选择的算法对其做 SVD,比较不同方法的效果(比如计算小气一直和对应左右奇异向量的 误差，效率等)，计算时间和所需存储量等，根据结果提出对算法的认识。 1.误差 在实验中，我们取 m=200，n=100，利用 orth()函数生成了正交矩阵、，再生成 了不同奇异值分布的奇异值矩阵，再通过 = ，计算出不同的待分解矩阵。 各矩阵奇异值分不如下表所示 序号 奇异值个数 奇异值分布 1 10 10 1 1 2 20 20 1 1 3 30 30 1 1 4 40 40 1 1 5 50 50 1 1 6 60 60 1 1 7 70 70 1 1 8 80 80 1 1 9 90 90 1 1 10 100 100 1 1 经过 matlab 计算，我们得到了两种算法对奇异值的估计误差表，如下所示 序号 2 100 1 100 ii i svd lansvd 1 1.2135e-29 1.2609e-29 2 1.0336e-28 1.3584e-28 3 3.9976e-28 2.7581e-28 4 7.9768e-28 7.2521e-28 5 1.5711e-27 1.5930e-27 6 2.7966e-27 2.5232e-27 7 3.5970e-27 4.6720e-27 8 6.2227e-27 5.8535e-27 9 7.2128e-27 8.7958e-27 10 1.0021e-26 9.8885e-27 奇异向量的估计误差，如下所示 序号 2 2 1 r iis i uu r 2 2 1 r iis i vv r svd lansvd svd lansvd 1 2 1.6 2 1.6 2 1.6 2.2 1.6 2.2 3 2.1333 2.1333 2.1333 2.1333 4 2 1.8 2 1.8 5 2.16 2.16 2.16 2.16 6 1.8 1.8667 1.8 1.8667 7 2.3429 2.2286 2.3429 2.2286 8 2.15 2.0275 2.15 2.2057 9 2 2.2075 2 2.2075 10 2 1.92 2 1.92 通过对比上面两个误差表，我们发现两种方法在对奇异向量矩阵的估计误差相近， 并没有明显的优异性，在估计奇异值上 svd 方法略强于 lansvd 方法。 他们的运行时间以及使用内存情况如下所示， svd 方法 lansvd 方法 由以上两数据对比可知，svd 方法效率明显优于 lansvd 方法，而二者的内存使用情 况相近。 综上， 我们不难得到如下结论， 两种方法在内存使用以及估计精度上没有明显的差 别，而在效率上 svd 方法远远优于 lansvd 方法。 三、图像压缩 与传统图像压缩的概念相同， SVD 方法也是舍弃部分图像数据从而达到压缩的效果， 这里我们舍弃的数据是奇异值较小奇异向量所包含的数据，具体方法如下所示: 1.将图像 A 做奇异值分解，即 0 00 AUV 2.将中较小的奇异值舍去，并且抛弃其所对应的奇异向量舍弃，并将剩余的数据 存储，形成压缩的图片信息 我们选择的 5 幅数字图像如下所示： (1)黑白网格图 (2)横向渐变灰度图 (3)斜向渐变灰度图 (4)Lena 灰度图 (5)白噪声灰度图 我们选择抛弃小于最大奇异值的 1%所对应的的奇异值与奇异向量如下所示： (1)黑白网格图 (2)横向渐变灰度图 (3)斜向渐变灰度图 (4)Lena 灰度图 (5)白噪声灰度图 通过对比以上五张数字图像， 直观上我们发现， 形状越复杂的图像压缩后的效果越 差，我们利用峰值信噪比来评价图片压缩后的质量 序号 保留特征向量数 2 2 1 255 10 log / N i PSNR OringinAfterN 1 3 639.1492 2 2 636.4492 3 4 682.0049 4 55 70.8480 5 657 58.7780 通过上表中的数据对比， 我们就可以发现越复杂的图像越难压缩， 且压缩后的质量 很差。 比较 SVD 方法与传统的 JPEG 方法可知， 再存储简单图像时 SVD 方法可能会获得更 高的压缩率，因为他不需要想 JPEG 记录图像上每块儿的信息，而是直接对整体图像进 行变换存储，而对于较复杂的图像时，由于 SVD 方法需要记录很多的奇异向量所以他 的压缩率可能会低于有固定正交向量矩阵的 JPEG 方法。 附： generator.m 为第二部分中生成计算矩阵的程序 my_svd.m 为利用不同房法求解矩阵奇异值分解的程序 fig_process. 为处理图像的程序 lansvd.m 为 Propack 中的奇异值分解程序