7.1 线性代数进行图像处理

一、图像的奇异值分解介绍

AAA 的奇异值定理就是 ATAA^TAATAAATAA^TAAT 的特征值定理。 这是本章的内容的预览。AAA 有两个奇异向量的集合(ATAA^TAATAAATAA^TAAT 的特征向量),有一个是正奇异值的集合(因为 ATAA^TAATAAATAA^TAAT 有相同的正特征值)。通常 AAA 是长方形的矩阵,但是 ATAA^TAATAAATAA^TAAT 是方阵,对称的且至少是半正定的。
奇异值分解(The Singular Value Decomposition:SVD)可以将任意的矩阵分解成简单的片(pieces)。 每个片都是一个列向量乘一个行向量,一个 m×nm\times nm×n 的矩阵有 m×nm\times nm×n 个元素(如果表示的是一个图像则会是很大的数字),但是列和行仅有 m+nm+nm+n 个分量,远少于 m×nm\times nm×n。这些列乘行 (colum)(row) 的小片是可以非常快速处理的全尺寸大小矩阵(full size matrices)—— 它们仅需要 m+nm+nm+n 个数。
不寻常的是,SVD 的图像处理应用是在它所依赖的代数矩阵之前出现的,下面会从只包含一个或两个片的简单图像开始。现在我们将图像想象成一个大的长方形矩阵,元素 aija_{ij}aij 表示的是图像所有像素的灰度值,可以将像素想象成一个小的正方形,从左下角开始向右 iii 步和向上 jjj 步,它的灰度值是一个数字(通常的范围是 0≤aij≤256=280\le a_{ij}\le256=2^80aij256=28),一个纯白的像素是 aij=255=11111111a_{ij}=255=11111111aij=255=11111111,计算机将 255255255 用二进制表示时,这个数字就是 888111
目前我们已经知道一个图像有 m×nm\times nm×n 个像素,每个像素使用 888 位(000111)表示它的灰度,变成了一个 m×nm\times nm×n 的矩阵,其中的每个元素都有 256256256 种的可能值。
简单的说,一个图像就是一个大的矩阵,为了完美复制它,我们需要 8(m)(n)8(m)(n)8(m)(n) 比特的信息。高分辨率的电视通常是 m=1080m=1080m=1080n=1920n=1920n=1920,每秒有 242424 帧,如果是彩色的,会有 333 个颜色度,那么需要每秒传输 (3)(8)(24)(1920)(1080)(3)(8)(24)(1920)(1080)(3)(8)(24)(1920)(1080) 位,这个成本太高了所以也不太会这么做,发射器跟不上要显示的速度。
如果压缩做的很好的话,那么我们是无法看出和原始图像的区别的。图像的边缘(灰度的突然变化)是很难压缩的部分。
如果每个 aija_{ij}aij 都是独立的随机数字,那么就不可能在压缩方面取得大的成果。我们完全依赖这样的一个事实:相邻的像素通常由相似的灰度,当跨越边缘时就会出现突然的跳变。卡通图片要比现实图片更容易压缩,这是因为它到处都是边缘。
对于视频来说,数字 aija_{ij}aij 不会在帧与帧之间变化太大,我们只传输小的变化。这就是 H.264 视频压缩标准的差分编码。我们使用线性代数(和非线性的 “量化quantization”,这是计算机中高效的转换成正数的一步)压缩每个变化矩阵。
我们每天看到的自然图像完全适合且可以进行压缩处理。

二、低秩图像(例子)

最容易压缩的图像是全黑或全白或全是一个常数灰度 ggg,矩阵 AAA 的每个元素的数都是同样的 gggaij=ga_{ij}=gaij=g,当 g=1g=1g=1m=n=6m=n=6m=n=6 时,下面是图像处理中关于 SVD 核心理论的一个极端的例子:

例1不是发送 A=[111111111111111111111111111111111111]而是发送 A=[111111][111111]\pmb{不是发送}\,A=\begin{bmatrix}1&1&1&1&1&1\\1&1&1&1&1&1\\1&1&1&1&1&1\\1&1&1&1&1&1\\1&1&1&1&1&1\\1&1&1&1&1&1\end{bmatrix}\kern 5pt\pmb{而是发送}\,A=\begin{bmatrix}1\\1\\1\\1\\1\\1\end{bmatrix}\begin{bmatrix}1&1&1&1&1&1\end{bmatrix}不是发送A=111111111111111111111111111111111111而是发送A=111111[111111]363636 个数字变成了 121212 个数字,如果是 300×300300\times300300×300 个像素,那么就是 900009000090000 个数变成 600600600 个数。如果我们提前定义好全一向量 x\boldsymbol xx,我们只需要发送一个数字,这个数字就是常数灰度 ggg,它乘上 xxT\boldsymbol x\boldsymbol x^TxxT 得到这个矩阵。
当然第一个例子很极端,但它阐明了一个重要观点,如果能够预先定义一些特殊的向量如全一向量 x=ones\boldsymbol x=\textrm{\pmb{ones}}x=ones,那么图像处理将会变得非常高效,本质上是预选基(preselected bases)(如傅里叶基通过 FFT 进行加速)和自适应基(adaptive bases) 之间的较量。SVD 通过图像本身生成基 —— 这个是自适应的,它的计算成本可能比较高。
并不是说 SVD 总是或者通常在实际上是最有效的算法,下面这些例子的目的只是为了介绍而不是实际生产使用。

例2“ace  三色旗”法国国旗  A意大利国旗 A德国国旗  AT不要发送 A=[aacceeaacceeaacceeaacceeaacceeaaccee]发送 A=[111111][aaccee]\begin{array}{l}“\pmb{\textrm{ace}\,\,三色旗}”\\法国国旗\,\,A\\意大利国旗\,A\\德国国旗\,\,A^T\end{array}\kern 5pt\pmb{不要发送}\,A=\begin{bmatrix}a&a&c&c&e&e\\a&a&c&c&e&e\\a&a&c&c&e&e\\a&a&c&c&e&e\\a&a&c&c&e&e\\a&a&c&c&e&e\end{bmatrix}\kern 5pt\pmb{发送}\,A=\begin{bmatrix}1\\1\\1\\1\\1\\1\end{bmatrix}\begin{bmatrix}a&a&c&c&e&e\end{bmatrix}ace三色旗法国国旗A意大利国旗A德国国旗AT不要发送A=aaaaaaaaaaaacccccccccccceeeeeeeeeeee发送A=111111[aaccee]这个虽说有 333 种颜色,但是它的秩仍然是 111,我们仍然可以用一列乘上一行。这 363636 个元素甚至可能都不相同,但是只要是秩一的模式 A=u1v1TA=\boldsymbol u_1\boldsymbol v_1^TA=u1v1T 都可以压缩。但是如果秩提升至 r=2r=2r=2,就需要 u1v1T+u2v2T\boldsymbol u_1\boldsymbol v_1^T+\boldsymbol u_2\boldsymbol v_2^Tu1v1T+u2v2T,例如:

例3内嵌的方块矩阵 A=[1011]\boxed{A=\begin{bmatrix}1&0\\1&1\end{bmatrix}}A=[1101] 等于 A=[11][11]−[10][01]\boxed{A=\begin{bmatrix}1\\1\end{bmatrix}\begin{bmatrix}1&1\end{bmatrix}-\begin{bmatrix}1\\0\end{bmatrix}\begin{bmatrix}0&1\end{bmatrix}}A=[11][11][10][01]

这里矩阵 AAA 中的 111000 可以为块状的 1′s1's1s0′s0's0s,它的秩仍为 222,我们同样只需要两项 u1v1T+u2v2T\boldsymbol u_1\boldsymbol v_1^T+\boldsymbol u_2\boldsymbol v_2^Tu1v1T+u2v2T。一个 6×66\times66×6 的图像可以被压缩成 242424 个数字,一个 N×NN\times NN×N 的图像(N2N^2N2 个数字)能够压缩成来自向量 u1,v1,u2,v2\boldsymbol u_1,\boldsymbol v_1,\boldsymbol u_2,\boldsymbol v_2u1,v1,u2,v2 这四个向量的 4N4N4N 个数字。
这里 u′s\boldsymbol u'susv′s\boldsymbol v'svs 是最佳的选择吗?这个并不是 SVD 的选择!注意到 u1=(1,1)\boldsymbol u_1=(1,1)u1=(1,1)u2=(1,0)\boldsymbol u_2=(1,0)u2=(1,0) 并不是正交的,v1=(1,1)\boldsymbol v_1=(1,1)v1=(1,1)v2=(0,1)\boldsymbol v_2=(0,1)v2=(0,1) 也不正交。理论表明:正交性能够使得第二项 c2u2v2Tc_2\boldsymbol u_2\boldsymbol v_2^Tc2u2v2T 更小。(SVD 按照重要性的顺序选择秩一项。)
如果 AAA 的秩远大于 222,就比如实际图像,那么 AAA 就会是很多秩一项相加,我们希望最小的项足够小 —— 它们可以被丢弃而不会影响视觉的质量,这时图像压缩就变为有损的了,但是好的图像压缩人的视觉几乎是无法察觉的。
问题变成了:SVD 的正交基要如何选择?

三、SVD 的特征向量

大部分图像的特征向量并不正交,而且特征向量 x1,x2\boldsymbol x_1,\boldsymbol x_2x1,x2 只提供一组向量,我们想要两组向量(u′s\boldsymbol u'susv′s\boldsymbol v'svs)。这两个问题的解决方法正是 SVD 的思想:
使用 AATAA^TAAT 的特征向量 u\boldsymbol uuATAA^TAATA 的特征向量 v\boldsymbol vv
因为 AATAA^TAATATAA^TAATA 自动就是对称的(通常并不相等!),所以特征向量 u′s\boldsymbol u'sus 是一组正交向量,特征向量 v′s\boldsymbol v'svs 是另一组正交向量,我们可以并且需要将它们单位化:∣∣ui∣∣=1||\boldsymbol u_i||=1∣∣ui∣∣=1∣∣vi∣∣=1||\boldsymbol v_i||=1∣∣vi∣∣=1,这样我们的秩 222 矩阵就是 A=σ1u1v1T+σ2u2v2TA=\sigma_1\boldsymbol u_1\boldsymbol v_1^T+\sigma_2\boldsymbol u_2\boldsymbol v_2^TA=σ1u1v1T+σ2u2v2T。其中 σ1\sigma_1σ1σ2\sigma_2σ2 这些数字的大小将决定它们在压缩过程中是否可以被忽略。我们保留大的 σ′s\sigma'sσs,丢弃小的 σ′s\sigma'sσs
SVD 中的 u′s\boldsymbol u'sus 称为左奇异值向量(left singular vectors),它们是 AATAA^TAAT 的单位特征向量;v′s\boldsymbol v'svs 称为右奇异值向量(right singular vectors),它们是 ATAA^TAATA 的单位特征向量。其中 σ′s\sigma'sσs 就是奇异值(singular values),等于 AATAA^TAATATAA^TAATA 特征值的平方根:

SVD 的选择 AATui=σi2uiATAvi=σi2viAvi=σiui(7.1.1){\color{blue}AA^T\boldsymbol u_i=\sigma_i^2\boldsymbol u_i\kern 10ptA^TA\boldsymbol v_i=\sigma_i^2\boldsymbol v_i\kern 10ptA\boldsymbol v_i=\sigma_i\boldsymbol u_i}\kern 45pt(7.1.1)AATui=σi2uiATAvi=σi2viAvi=σiui(7.1.1)

在例 3 中(嵌入的方块矩阵),下面是对称矩阵 AATAA^TAATATAA^TAATAAAT=[1011][1101]=[1112]ATA=[1101][1011]=[2111]AA^T=\begin{bmatrix}1&0\\1&1\end{bmatrix}\begin{bmatrix}1&1\\0&1\end{bmatrix}=\begin{bmatrix}\pmb1&\pmb1\\\pmb1&\pmb2\end{bmatrix}\kern 20ptA^TA=\begin{bmatrix}1&1\\0&1\end{bmatrix}\begin{bmatrix}1&0\\1&1\end{bmatrix}=\begin{bmatrix}\pmb2&\pmb1\\\pmb1&\pmb1\end{bmatrix}AAT=[1101][1011]=[1112]ATA=[1011][1101]=[2111]它们的特征值都是 111,所以 λ1λ2=1\lambda_1\lambda_2=1λ1λ2=1,它们的迹(对角元素之和)都是 333det⁡[1−λ112−λ]=λ2−3λ+1=0得到λ1=3+52和λ2=3−52\det\begin{bmatrix}1-\lambda&1\\1&2-\lambda\end{bmatrix}=\lambda^2-3\lambda+1=0\kern 5pt得到\kern 5pt\lambda_1=\frac{3+\sqrt5}{2}\kern 5pt和\kern 5pt\lambda_2=\frac{3-\sqrt5}{2}det[1λ112λ]=λ23λ+1=0得到λ1=23+5λ2=235λ1\lambda_1λ1λ2\lambda_2λ2 的平方根是 σ1=5+12\sigma_1=\displaystyle\frac{\sqrt5+1}{2}σ1=25+1σ2=5−12\sigma_2=\displaystyle\frac{\sqrt5-1}{2}σ2=251σ1σ2=1\sigma_1\sigma_2=1σ1σ2=1。最接近 AAA 的秩一矩阵是 σ1u1v1T\sigma_1\boldsymbol u_1\boldsymbol v_1^Tσ1u1v1T,误差仅为 σ2≈0.6\sigma_2\approx0.6σ20.6,这个是理论最小值。
AATAA^TAATATAA^TAATA 标准正交特征向量是u1=[1σ1]u2=[σ1−1]v1=[σ11]v2=[1−σ1]全部除以1+σ12(7.1.2)\boldsymbol u_1=\begin{bmatrix}1\\\sigma_1\end{bmatrix}\kern 10pt\boldsymbol u_2=\begin{bmatrix}\sigma_1\\-1\end{bmatrix}\kern 10pt\boldsymbol v_1=\begin{bmatrix}\sigma_1\\1\end{bmatrix}\kern 10pt\boldsymbol v_2=\begin{bmatrix}1\\-\sigma_1\end{bmatrix}\kern 10pt全部除以\kern 3pt\sqrt{1+\sigma_1^2}\kern 20pt(7.1.2)u1=[1σ1]u2=[σ11]v1=[σ11]v2=[1σ1]全部除以1+σ12(7.1.2)实际生活中这些计算都是由计算机完成的!(MATLAB 中可以使用 svd(A) 进行奇异值分解。)我们可验证一下 σ1u1v1T+σ2u2v2T\sigma_1\boldsymbol u_1\boldsymbol v_1^T+\sigma_2\boldsymbol u_2\boldsymbol v_2^Tσ1u1v1T+σ2u2v2T 是否可以正确的重构矩阵 AAA

A=[u1u2][σ1σ2][v1Tv2T]\boxed{A=\begin{bmatrix}\boldsymbol u_1&\boldsymbol u_2\end{bmatrix}\begin{bmatrix}\sigma_1&\\&\sigma_2\end{bmatrix}\begin{bmatrix}\boldsymbol v_1^T\\\boldsymbol v_2^T\end{bmatrix}}A=[u1u2][σ1σ2][v1Tv2T] 或更简单的 A[v1v2]=[σ1u1σ2u2])(7.1.3)\boxed{A\begin{bmatrix}\boldsymbol v_1&\boldsymbol v_2\end{bmatrix}=\begin{bmatrix}\sigma_1\boldsymbol u_1&\sigma_2\boldsymbol u_2\end{bmatrix})}\kern 30pt(7.1.3)A[v1v2]=[σ1u1σ2u2])(7.1.3)

重点:关键点是图像并不都是趋于低秩,大部分图像都是满秩的,但是它们有低有效秩(low effective rank)。这意味着:很多奇异值都很小,它们可以设置成零,我们传输的是一个低秩的近似。

例4】假设旗帜是两种不同颜色的三角形,左下角的三角形都是 111,右上角的三角形都是 000,主对角线都是 111。下面是当 n=4n=4n=4 时的图像矩阵,它是满秩的 r=4r=4r=4,所以它可逆:三角形旗帜矩阵A=[1000110011101111]和A−1=[1000−11000−11000−11]\begin{array}{l}\pmb{三角形}\\\pmb{旗帜矩阵}\end{array}\kern 10ptA=\begin{bmatrix}\pmb1&0&0&0\\\pmb1&\pmb1&0&0\\\pmb1&\pmb1&\pmb1&0\\\pmb1&\pmb1&\pmb1&\pmb1\end{bmatrix}\kern 5pt和\kern 5ptA^{-1}=\begin{bmatrix}\kern 7pt1&\kern 7pt0&\kern 7pt0&0\\-1&\kern 7pt1&\kern 7pt0&0\\\kern 7pt0&-1&\kern 7pt1&0\\\kern 7pt0&\kern 7pt0&-1&1\end{bmatrix}三角形旗帜矩阵A=1111011100110001A1=1100011000110001由于矩阵满秩,所以 AAA 有一组 nnn 个奇异值 σ\sigmaσ(都是正数),SVD 会产生 nnn 个秩一项 σiuiviT\sigma_i\boldsymbol u_i\boldsymbol v_i^TσiuiviT,完美的重构需要全部的 nnn 项。
在压缩过程中小的 σ\sigmaσ 可以被丢弃,它们不会明显影响图片的质量。我们要理解并画出当 n=4n=4n=4 时的这些 σ\sigmaσ 以及更大的 nnn。注意例 3 是例 4 这个三角形矩阵当 n=2n=2n=2 时的特例。
下面我们手动计算,从 AATAA^TAAT 开始(计算机处理会不同):AAT=[1111122212331234]和(AAT)−1=(A−1)TA−1=[2−100−12−100−12−100−11](7.1.4)AA^T=\begin{bmatrix}1&1&1&1\\1&2&2&2\\1&2&3&3\\1&2&3&4\end{bmatrix}\kern 10pt和\kern 10pt(AA^T)^{-1}=(A^{-1})^TA^{-1}=\begin{bmatrix}\kern 7pt\pmb2&-1&\kern 7pt0&\kern 7pt0\\-1&\kern 7pt\pmb2&-1&\kern 7pt0\\\kern 7pt0&-1&\kern 7pt\pmb2&-1\\\kern 7pt0&\kern 7pt0&-1&\kern 7pt\pmb1\end{bmatrix}\kern 15pt(7.1.4)AAT=1111122212331234(AAT)1=(A1)TA1=2100121001210011(7.1.4)这个 −1,2,−1-1,2,-11,2,1 的逆矩阵被引进是因为它的特征值的形式是 2−2cos⁡θ2-2\cos\theta22cosθ,所以我们可以得到 AATAA^TAATλ′s\lambda'sλsAAAσ′s\sigma'sσsλ=12−2cos⁡θ=14sin⁡2(θ/2)得到σ=λ=12sin⁡(θ/2)(7.1.5)\lambda=\frac{1}{2-2\cos\theta}=\frac{1}{4\sin^2(\theta/2)}\kern 10pt得到\kern 10pt\sigma=\sqrt\lambda=\frac{1}{2\sin(\theta/2)}\kern 25pt(7.1.5)λ=22cosθ1=4sin2(θ/2)1得到σ=λ=2sin(θ/2)1(7.1.5)nnn 个不同的角度 θ\thetaθ000 ~ πππ 之间等距离的,注意这个三对角矩阵最后一个元素是 111 而不是 222,此时等距的点有 2n2n2n 个,若为 222 时等距点有 nnn 个,这个例子非常特殊:θ=π2n+1,3π2n+1,⋯ ,(2n−1)π2n+1(n=4 时包括 θ=3π9,此时 2sin⁡θ2=1)\theta=\frac{π}{2n+1},\frac{3π}{2n+1},\cdots,\frac{(2n-1)π}{2n+1}\kern 10pt\Big(n=4\,时包括\,\theta=\frac{3π}{9},此时\,2\sin\frac{\theta}{2}=1\Big)θ=2n+1π,2n+13π,,2n+1(2n1)π(n=4时包括θ=93π,此时2sin2θ=1)这个特殊案例得到了当 n=4n=4n=4AATAA^TAAT 的一个特征值 λ=1\lambda=1λ=1,所以 σ=λ=1\sigma=\sqrt\lambda=1σ=λ=1AAA 的一个奇异值,当向量 u=(1,1,0,−1)\boldsymbol u=(1,1,0,-1)u=(1,1,0,1) 时,有 AATu=uAA^T\boldsymbol u=\boldsymbol uAATu=u,这是一个真正的特例。
重点是画出 AAAnnn 个奇异值的图形,这些数字是逐渐减小的(不像 AAA 的特征值,都是 111),但是减小的不是很陡峭。所以 SVD 只是适度的压缩这个三角形旗帜。对于希尔伯特矩阵压缩效果非常好。

在这里插入图片描述

四、例题

例5】MATLAB 指令 A = rand(20, 40) 和 B = randn(20, 40) 会生成 20×40 的随机矩阵:AAA 中的元素是 000111 之间均匀分布的;B 中的元素是 “钟形” 正态分布。使用 MATLAB 中的 svd 命令,画出奇异值 σ1σ_1σ1σ20σ_{20}σ20
解: MATLAB 代码如下:

% 生成矩阵并计算奇异值
A = rand(20, 40);
s = svd(A);

% 绘制半对数坐标图
figure;
semilogy(s, '-o', 'LineWidth', 1.5, 'MarkerSize', 6); % 这个函数是绘制半对数坐标的函数,横轴是线性,纵轴是对数;
													  %-o 是圆圈标记;MarkerSize 是这个圆圈标记的大小
xlabel('序号 (k)');
ylabel('\sigma_k (对数刻度)');
title('A=rand(20,40)的奇异值分布');
axis tight;											  % 自动调整坐标轴范围,使数据内容紧密贴合坐标轴边界
grid on;

% 生成20×40的随机矩阵(元素服从标准正态分布)
B = randn(20, 40);

% 计算奇异值分解(Singular Value Decomposition)
[U, S, V] = svd(B);

% 提取奇异值(S矩阵的对角线元素)
singular_values = diag(S);

% 绘制奇异值分布图(对数坐标系)
figure;
semilogy(singular_values, 'bo-', 'LineWidth', 1.5);	% bo 表示蓝色的圆圈,- 代表实线
grid on;
xlabel('序号 (k)');
ylabel('\sigma_k(对数刻度)');
title('B=randn(20,40)的奇异值分布');
axis tight;

运行结果如下:

图片1地址![](https://img-blog.csdnimg.cn/图片2地址
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值