基于数字实例的奇异值分解教学研究

基于数字实例的奇异值分解教学研究
刘智慧;付丽华
【摘要】考虑到奇异值分解涉及的理论性较强,试图从奇异值分解的基本原理出发,通过图像压缩与恢复和图像降噪2个具体的实例,结合MATLAB数学软件来阐述奇异值分解的作用,注重应用奇异值分解理论解决相关的实际问题,以期加深学生对奇
异分解理论的理解,领会奇异值分解的意义和价值.
【期刊名称】《海南大学学报（自然科学版）》
【年(卷),期】2015(033)004
【总页数】5页(P391-395)
【关键词】奇异值分解;图像压缩恢复;图像降噪
【作者】刘智慧;付丽华
【作者单位】中国地质大学数理学院,武汉,湖北430074;中国地质大学数理学院,武汉,湖北430074
【正文语种】中文
【中图分类】O151
矩阵的奇异值分解是最重要的矩阵分解之一，具有较高的理论价值和极其广泛的实际应用[1-11]，其理论由Beltrami和Jordan分别于1873年和1874年独立创建. 由于奇异值分解理论较为抽象深奥，学生在学习时比较容易陷入枯燥的理论证明而失去对其学习的热情和兴趣. 在大学课程《线性代数》、《计算方法》和《矩阵论》等课程中也很少介绍奇异值分解有关的应用背景和实际应用. 如果教师在讲授奇异
值分解定理时，通过一些实际生活中的具体例子来说明矩阵奇异值分解的作用，将在一定程度上调动学生的学习积极性，激发其学习兴趣，从而使其较好地掌握矩阵的奇异值分解.通过图像压缩与恢复和图像降噪2个具体实例，结合MATLAB数学软件说明奇异值分解的作用，以求加深学生对奇异值分解的真正理解，拓展奇异值分解理论的进一步应用.
定理1 奇异值分解定理[1] 设矩阵A则存在m阶酉矩阵U=[u1，
u2,…,ur,ur+1,…,um]和n阶酉矩阵V=[v1,v2,…,vr,vr+1,…,vn]，使得
其中，(·)H表示复共轭转置运算，∑=diag{σ1,σ2,…,σr}且σ1≥σ2≥…≥σr>0，称σi为A的奇异值，列向量ui(vi)分别称为A对应奇异值σi的左(右)奇异向量.
从定理1中可以得到，任意一个比较复杂的矩阵A总可以用阶数更小更简单的一些子矩阵ui和vi相乘来表示.
在实际应用中，得到的数据矩阵往往会带有噪声或扰动，这样就会使得原本是低秩矩阵的奇异值全部为非零. 为此，可对所得数据矩阵进行奇异值分解，利用其矩阵前若干个较大奇异特征值和相对应的左右奇异向量来重构原数据矩阵.
定理2[2-3] 设矩阵A的奇异值分解为则为矩阵A在Frobenius范数意义下的最佳近似矩阵，即
从定理2中可以得到，在Frobenius范数意义下，矩阵A的所有秩为k的近似矩阵中，Ak最接近原始矩阵A，可以理解为Ak对矩阵A的重构，即Ak≈A.
通过一个数据矩阵的奇异值分解来考察奇异值对矩阵的影响.
例1 求矩阵的奇异值分解.
利用MATLAB命令求出奇异值分解的酉矩阵和对角矩阵，相应的MATLAB程序为：
A=[1 2 3;5 1 2;7 1 3];
[U,D,V]=svd(A);
运行结果如下：
U =
-0.264 3 0.964 3 0.016 7
-0.560 0 -0.139 4 -0.816 6
-0.785 2 -0.225 2 0.576 9
D =
9.747 7 0 0
0 2.813 8 0
0 0 0.255 2
V =
-0.878 2 - 0.465 2 -0.111 0
-0.192 2 0.555 8 -0.808 8
-0.437 9 0.688 9 0.577 6
故A的奇异值分解为
A= UDVH=
.
从矩阵A的奇异值分解可以看到：对角矩阵∑中的第三个奇异值σ3=0.255 2相比于另外2个奇异值较小，将其忽略，即令σ3=0. 计算前2个大的奇异值σ1，σ2和相对应的左右奇异向量u1,u2,v1,v2构成的矩阵A2，即
比较矩阵A2和A，有A2≈A.从例1可以得到，通过奇异值分解，利用原数据矩阵的前2个较大特征值和相应的左右奇异向量可以重构原数据矩阵A.
2.1 图像压缩与恢复现有一幅卫星图像，假设此图像含n×n个像素，即这张图像包含n2个数据，现将此图像发送到地球，如果将n2个数据一起传送，那么数据量往往很大. 因此，希望能否只传送n2个数据中的较少数据，同时又能在地球接
收端可以通过接收的数据重构原图像？回答是肯定的.
将此卫星图像看作n×n的矩阵A，假设已知矩阵A的奇异值分解，即设前
k(k<n)个奇异值σi(i=1,…,k)很大，其他n-k个奇异值很小. 那么，选择前k个大
的奇异值σi(i=1,…,k)及矩阵U和V的前k列发送给地球上的接收端，即矩阵A近似表示成
=Ak ，
式(1)表明，可以用k(2n+1)个数据代替原来的n2个图像数据，即用矩阵Ak对应的图像恢复原始图像A. 进一步，定义图像的压缩比为n2：k(2n+1)，当
k(2n+1)<n2时，可以达到图像压缩的目的[4].
用具体例子说明基于奇异值分解的图像压缩与恢复. 图1(a)表示含128×128像素
的Lena图像，将此图像可以看作是一个128×128阶的矩阵A. 对此图像矩阵A
进行奇异值分解可得一条奇异值曲线，如图2所示. 从图2可以看到，随着k的增加，奇异值σk迅速减小.k=1时，奇异值σ1=13 257，当k≥20时，奇异值
σk≤378.
图1表示利用奇异值分解对Lena图像的压缩与恢复. 图1(b)，(c)，(d)分别表示矩阵A20,A40,A60所对应的图像，其中A20,A40,A60分别表示A矩阵前20个，
40个，60个大的奇异值及与这些奇异值相对应的左(右)奇异向量构成的矩阵.
从图1可知，利用奇异值分解对图像进行压缩与恢复时，k的取值越大，恢复的图像效果越好，但压缩率会越小；k的取值越小，压缩率会越大，但恢复的图像效果越不理想. 例如，在图1(b)中，k=20时，图像的压缩率为3.19∶1，压缩率较大，但恢复的图像与原始图像1(a)相差较大，图像的恢复效果不理想. 在图3(c)中，
k=40时，恢复的图像与原始图像1(a)非常接近，图像的压缩率为1.59∶1. 在图
3(d)中， k=60时，恢复的图像效果较好，与原始图像1(a)非常接近，但图像的压缩率只有1.06∶1. 所以，对图像1(a)而言，取k=40时，图像既有较高的压缩比
同时图像又得到了较好的恢复.
以上例子说明，利用矩阵的奇异值分解，选择矩阵前面部分较大奇异值而忽略较小奇异值可以实现图像的压缩与恢复.
图1相应的MATLAB程序为：
A= imread(Lena.bmp); % 读取128x128的Lena图像
A = double(A); %数据格式转换
[U,D,V]=svd(A); % 对数据矩阵进行奇异值分解
recv1=U(:,1:20)*D(1:20,1:20)*V(:,1:20); %提取前20个大的奇异值和奇异向量recv2=U(:,1:40)*D(1:40,1:40)*V(:,1:40); %提取前40个大的奇异值和奇异向量recv3=U(:,1:60)*D(1:60,1:60)*V(:,1:60); %提取前60个大的奇异值和奇异向量subplot(221);
imshow(A);
xlabel((a) 原始图像);
subplot(222);
imshow(mat2gray(recv1));
xlabel((b) k=20);
subplot(223);
imshow(mat2gray(recv2));
xlabel((c) k=40);
subplot(224);
imshow(mat2gray(recv3));
xlabel((d) k=60);
图2相应的MATLAB程序为：
A= imread(Lena.bmp); % 读取128x128的Lena图像
A = double(A); %数据格式转换
[U,D,V]=svd(A); % 对数据矩阵进行奇异值分解
for i=1:128 %提取奇异值
d(i)=D(i,i);
end
i=1:128;
plot(i,d(i),*r);
xlabel(k);
ylabel(\sigma_k);
2.2 图像降噪利用奇异值分解对图像进行降噪. 设有一张含15×11个像素的黑白
图像,见图3(a). 用15×11的矩阵A表示该图像，对于黑白图像，像素是白的就取数1，像素是黑色的就取数0. 现将此黑白图像添加均值为零，标准差为0.1的高
斯白噪声，得到一幅含噪声的图像，见图3(b)(用矩阵B表示).
为了对含噪声图像3(b)进行降噪，计算图像3(b)的奇异值(即计算矩阵B的奇异值)，其奇异值分别为：σ1=6.707 3，σ2=0.607 4，σ3=0.554 1，σ4=0.488 0，σ5=0.391 7，σ6=0.370 4，σ7=0.289 2，σ8=0.234 7，σ9=0.216 7，
σ10=0.110 1，σ11=0.084 9. 显然，第1个奇异值要比后面的奇异值大很多. 取
第一个奇异值对应的矩阵其中矩阵B1对应图3(c). 比较图3(a)，(b)和(c)可知，利用奇异值分解，选取矩阵较大奇异值，而舍去较小奇异值，可以对图像进行降噪，而且降噪效果较好.
图3相应的MATLAB程序为：
M=15; N=11;
a=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
b=[0 0 0 1 1 1 1 1 1 1 1 1 0 0 0];
A=[a a a b b b b b a a a];
e=normrnd(0,0.1,M,N);
B=A+e;
subplot(131);
imshow(A);
xlabel((a)原始图像);
subplot(132);
imshow(B);
xlabel((b)加噪图像);
[U,D,V]=svd(B);
C=U(:,1)*D(1,1)*V(:,1);
subplot(133);
imshow(C);
xlabel((c)降噪图像);
运用MATLAB数学软件将奇异值分解理论应用到解决图像压缩与恢复和图像降噪等实际问题的教学方式，使得原本枯燥难懂的抽象理论变得易于理解记忆，使得教学内容得到进一步的扩展. 同时，有利于培养学生学习数学的兴趣和热情，有利于培养其运用数学解决实际问题的能力. 此教学方式将对提高数学教学质量和培养复合型人才起到一定的积极作用.
【相关文献】
[1] 张贤达. 矩阵分析与应用[M]. 北京：清华大学出版社, 2004.
[2] 聂守平, 魏晓燕. 数字图像的奇异值分解[J]. 南京师大学报(自然科学版), 2001, 24(1): 9-61.
[3] 尹芳黎, 杨雁莹, 王传栋, 等. 矩阵奇异值分解及其在高维数据处理中的应用[J]. 数学的实践与认
识, 2011, 41(15): 171-176.
[4] 罗小桂. 矩阵奇异值分解(SVD)的应用[J]. 井冈山医专学报, 2005, 12(4): 133-135.
[5] 申卯兴, 郑武团. 矩阵的奇异值分解的应用[J]. 工科数学, 1996, 12(3)： 56-61.
[6] 王宏兴, 霍玉洪. 奇异值分解的教与学[J]. 淮南师范学院学报, 2015, 17(2): 125-127.
[7] 张磊, 彭伟才, 原春晖, 等. 奇异值分解降噪的改进方法[J]. 中国舰船研究, 2012, 7(5): 83-88.
[8] 刘波, 杨华, 张志强. 基于奇异值分解的图像去噪[J]. 微电子学与计算机, 2007, 24(11): 169-171.
[9] 黄君, 包晓敏, 杜晶晶. 基于HVS和奇异值分解的图像质量评价[J]. 工业控制计算机, 2012, 25(3): 39-40.
[10] 高媛赟, 温小荣, 林国忠, 等. 基于奇异值分解的中国资源一号O2C卫星数据的融合评价及应用[J]. 南京林业大学学报(自然科学版), 2015, 39(3): 29-33.
[11] 侯翔, 杨辉. 基于奇异值分解的遥感图像鲁棒水印算法探讨[J]. 测绘科学, 2015, 40(7): 123-126.。