清华大学高等数值分析(李津)第二次实践作业

奇异值分解算法分析及利用SVD 进行图像压缩

报告提纲：

一、 奇异值分解 二、 SVD 算法 三、 实验过程 四、 实验结果分析

五、 SVD 法图像压缩及结果分析

奇异值分解

定义设m n A R ?∈，T A A 的特征值的非负平方根称作A 的奇异值；A 的奇异值的全

体记作()A σ。

当A 为复矩阵m n C ?时，只需将T A A 改为H A A ，定义仍然成立。 定理（奇异值分解定理）设m n A R ?∈，则必存在正交矩阵

[]1,...,m m m U u u R ?=∈和[]1,...,n n n V v v R ?=∈

使得

00

0r

T r n r

r

U AV m r

-∑

??=??-?? 其中11(,...,),...0r r r diag σσσσ∑=≥≥>[2]。

当A 为复矩阵m n C ?时，只需将定理中,U V 改为酉矩阵，其它不变，定理仍然成立。

奇异值分解通常简称为SVD ，i σ是A 的奇异值，向量i u 和i v 分别是第i 个左奇异向量和第i 个右奇异向量。

SVD 算法

1. QR 迭代算法

（1）输入()m n A R m n ?∈≥及允许误差ε.

（2）二对角化：计算Householder 变换1,,,n P P ???12,,n H H -???使得

112()(),0T

n n B n

P P A H H m n -????????=??

-??

其中12

00n n B δγγδ????

?

?=?????

? ； 12:,n U PP P =???122:.n V H H H -=???

具体算法，书面作业中题目具体解释。

（3）收敛性检验： （i ）将所有满足

()

1j j j γεδδ-≤+

的j 置零；

（ii ）如果0,2,,j j n γ==???，则输出有关信息结束；否则，1:0γ=，确定正整数p q <，使得

10p q n γγγ+==???==，0j γ≠，p j q <≤；

（iii ）如果存在i 满足1p i q ≤≤-使得

i B δε∞≤，

则111:0,:,:,:0,:1i i i i x y l δγδγ+++=====，转步（iv ），否则转步（4）. （iv ）确定cos ,sin c s θθ==和σ使

0,c s x s c y σ??????

=??????-??????

//这也相对于,0T

c s y s c x σ??????

=????

??-??????

所以可以直接调用givens 变换算法得到

:,i l δσ+=:(,,);T U UG i i l θ=+

//这相当于

(1:;,)(1:;,)(1:;,);T

c s c s U n i i l U n i i l U n i i l s c s c -????

+=+=+????-????

（v ）如果l q i <-，则

1:i l x s γ++=，11:i l i l c γγ++++=，1:i l y δ++=，:1l l =+， 转步（iv ），否则转步（i ）.

（4）SVD 迭代：将书上的追赶法作用于二对角阵

1

12100p p p p q q B δγδγγδ+++??????

?

?=??

?

????? ， 得

11:T B P B Q =，

:(,,)p n p q U Udiag I P I --=, :(,,)p n p q V Vdiag I Q I --=

然后转步（3）。

2. 分而治之

分而治之算法的第一步也和QR 迭代算法一样，先将矩阵二对角化得到二对角矩阵B ；然后将求B 的SVD 问题分为两个子问题。 先将B 写成如下的形式：

1

1200k k

k B B e e B αβ?? ?= ? ??

?

其中(1)1k k B R -?∈，()()2n k n k B R -?-∈为上二对角矩阵，j e 是相应维数的向量中的第j 个单位向量，并且一般我们取/2k n =????。 现在假如我们已经求得了1B ，2B 的SVD 如下：

()11110T B Q D W =，2222T B Q D W =

并且令()11T l λ为1W 的最后一行，2T f 为2W 的第一行。那么将这些带入B ，我们

可以得到：

111

11

2222000

000100000

T

T

T k k k Q D W B l f W Q D ααλβ?????? ???= ? ????? ????

???

可以看出中间的矩阵形式上很简单，除了对角元素和第k 行上的元素，其它元素都为0；这里先把它的SVD 的计算问题放到最后，而假定已知其SVD 为：T S G ∑，将它带入上式，就可以得到B 的SVD 分解式：

T B Q W =∑

其中

1

20

00

1000Q Q S Q ??

?

= ? ??

?；1200W W G W ??= ??? 而计算1B ，2B 的SVD 时也可递归地应用这种办法，至到子问题足够得小。

3. 单边雅克比迭代法

虽然以上两种方法都可以对二对角矩阵求出具有很高的相对精度的奇异值和奇异值分解式；但是将一个稠密矩阵二对角化这个过程却可能会造成很大的相对误差[17]。所以所有基于先将矩阵二对角化的SVD 方法都不可能具有可靠的良好的相对精确度。以下所要介绍的Jacobi 方法和这些方法完全不同，它通过一系列平面旋转（一般文献中在介绍Jacobi 方法时为了方便也将平面旋转称为Jacobi 变换）；最终使得矩阵收敛于一个对角矩阵：

(1)()(0),0,1,2,...()k k k G G J k G G +==≡

其中每一个k J 设计成，使得对于一对选定的指标(,)k k p q ，有如下式子成立：

(1)(1),(())0k k k T k p q G G ++=

最终收敛时，令()0lim ,k k k k G G V J ∞∞=→∞

==∏，有G GV ∞=，且G ∞具有互相正交的列

向量，可以将G ∞写成?G U

∞=∑，其中∑为对角矩阵，?U 为正交矩阵；从而得到原矩阵的SVD 分解式。

在Jacobi 算法中，令000G B D =是最初的矩阵，m m m G B D =是1m G -经过一次Jacobi 旋转过得到的矩阵，其中m D 为对角矩阵，m B 是列向量范数全为1的矩阵。可以知道对G BD =进行单边的Jacobi 变换相当于对矩阵

T T H G G DB BD DAD ===进行双边的Jacobi 变换。

实验过程

针对以上三种SVD 算法用matlab 进行数值模拟，三个算法源自NAG matlab 工具箱（第24版，64bit ），分别为f08kb （QR 迭代），f08kd （分而治之）和f08kj （单边雅克比迭代）。

实验中选出10个矩阵（大小1200×960，将图像压缩应用的5张图片每张图片一分为二得到10个矩阵）对每种算法的内存占用量、运算时间、最小特征值误差，以及左特征向量和右特征向量的计算误差进行了计算，并得到了相应的曲线图。特征向量的误差通过()T U U I -和()T V V I -列向量范数的均方误差来评估。至于最小特征值误差，默认matlab 自带的svd 算法得到的结果是“准确”的，作为参考标准。事实上，无法得到准确的特征值，用matlab 自带svd 作为标准严格来说也是不准确的。

图片文件为附件中的image1.jpg ，image2.jpg ，image3.jpg ，image,4.jpg ，image,5.jpg 。

实验结果分析

下面对算法的内存使用，计算速度，最小特征值精度，奇异向量误差进行了比较，各图中红色曲线为QR 迭代法，绿色代表分而治之方法，蓝色代表单边雅克比方法。

1. 内存使用

很明显，分而治之内存占用最少，大概而QR 迭代的一半，雅克比最大，大概比QR 迭代多占用一半的内存。

2.计算时间

计算时间上，分而治之最快，QR迭代次之，单边雅克比法最慢，花费时间大概是传统QR的十倍，是分而治之的20倍。

3.最小特征值精度

从图中可以看到，传统QR算法和分而治之的最小奇异值几乎是完全一样的，而单边雅克比则不同，计算精度上没有明显的差别，但雅克比计算结果温度，而另外两种则抖动相对较大。对各误差求和得到三种算法的误差和分别为2.8363e-15，3.2270e-15和1.9147e-15，相比之下，雅克比算法精度最高，传统QR算法次之，分而治之精度最差，原因在

报告第一部分已经提及，但三者差别并不明显。由于计算过程中的系统

精度为 2.2204e-16，显然实验得到的最小奇异值的误差大概在系统精度

的5~10倍。

4.左奇异向量误差

奇异向量误差上来看，分而治之最小，这与分而治之的算法原理一致，矩阵越

小计算的正交阵正交度越高。而雅克比的奇异向量的误差相对来说就要大得多，大概为传统QR迭代的5倍。

5.右奇异向量误差

右奇异向量误差的规律和左奇异向量一致，但总体误差更小，因为右奇异向量就是计算中的特征向量，且矩阵阶数小，而左奇异向量是通

过奇异值和特征向量反推得到的，故精度要差。

总体来说，传统QR方法程序简单，但占用内存大，在计算精度和运算速度方面居中。

分而治之的方法在内存占用和运算速度方面占有绝对优势，计算精度方面，奇异值精度没有传统QR好，但左右奇异向量的正交性要好于

传统QR方法。

单边雅克比方法的奇异值计算精度方面有优势，内存占用居中，但奇异向量的正交性差，运算速度相比前两种方法慢一个数量级。

实际应用中，根据需要选择特定的SVD算法。

注：图像数据在附件中。

SVD法图像压缩

实验中所用照片为附件中image1.jpg，image2.jpg，image3.jpg，image,4.jpg，image,5.jpg。首先将彩色图片灰度化，得到两维的矩阵（1200×1920），采用matlab 自带的svd算法进行奇异值分解。取前l项奇异值和奇异向量的乘积求和得到的图像成为前l项压缩，l取20、80、300，给出压缩后的图像，并进行比较。以狗的照片为例分析，另外4张图像类似，见附件。

可以看到，l=20时有明显的不清晰，就连光滑的背景都有明显的杂点；l=80时已经很清晰了，但放大看仍有噪点和部分细节分辨率下降；而l=300时就几乎和原图的清晰度一样了。

SVD压缩事实上是一种有损压缩，压缩过程中丢掉了图片的大量信息，这些信息有可能是一些重要的细节。

附件清单

1.程序清单（main.m，SVDfun.m）

2.原彩色图片五张

3.实验得到的图片20张（每个图4张）

4.计算部分中的图像数据（matlab data.mat）