一、图像的奇异值分解介绍
AAA 的奇异值定理就是 ATAA^TAATA 和 AATAA^TAAT 的特征值定理。 这是本章的内容的预览。AAA 有两个奇异向量的集合(ATAA^TAATA 和 AATAA^TAAT 的特征向量),有一个是正奇异值的集合(因为 ATAA^TAATA 和 AATAA^TAAT 有相同的正特征值)。通常 AAA 是长方形的矩阵,但是 ATAA^TAATA 和 AATAA^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^80≤aij≤256=28),一个纯白的像素是 aij=255=11111111a_{ij}=255=11111111aij=255=11111111,计算机将 255255255 用二进制表示时,这个数字就是 888 个 111。
目前我们已经知道一个图像有 m×nm\times nm×n 个像素,每个像素使用 888 位(000 或 111)表示它的灰度,变成了一个 m×nm\times nm×n 的矩阵,其中的每个元素都有 256256256 种的可能值。
简单的说,一个图像就是一个大的矩阵,为了完美复制它,我们需要 8(m)(n)8(m)(n)8(m)(n) 比特的信息。高分辨率的电视通常是 m=1080m=1080m=1080 且 n=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 的每个元素的数都是同样的 ggg:aij=ga_{ij}=gaij=g,当 g=1g=1g=1 且 m=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 中的 111 和 000 可以为块状的 1′s1's1′s 和 0′s0's0′s,它的秩仍为 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'su′s 和 v′s\boldsymbol v'sv′s 是最佳的选择吗?这个并不是 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'su′s 和 v′s\boldsymbol v'sv′s)。这两个问题的解决方法正是 SVD 的思想:
使用 AATAA^TAAT 的特征向量 u\boldsymbol uu 和 ATAA^TAATA 的特征向量 v\boldsymbol vv。
因为 AATAA^TAAT 和 ATAA^TAATA 自动就是对称的(通常并不相等!),所以特征向量 u′s\boldsymbol u'su′s 是一组正交向量,特征向量 v′s\boldsymbol v'sv′s 是另一组正交向量,我们可以并且需要将它们单位化:∣∣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'su′s 称为左奇异值向量(left singular vectors),它们是 AATAA^TAAT 的单位特征向量;v′s\boldsymbol v'sv′s 称为右奇异值向量(right singular vectors),它们是 ATAA^TAATA 的单位特征向量。其中 σ′s\sigma'sσ′s 就是奇异值(singular values),等于 AATAA^TAAT 和 ATAA^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^TAAT 和 ATAA^TAATA:AAT=[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,它们的迹(对角元素之和)都是 333:det[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−λ]=λ2−3λ+1=0得到λ1=23+5和λ2=23−5λ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=25−1 且 σ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σ2≈0.6,这个是理论最小值。
AATAA^TAAT 和 ATAA^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=[σ1−1]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=1111011100110001和A−1=1−10001−10001−10001由于矩阵满秩,所以 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=(A−1)TA−1=2−100−12−100−12−100−11(7.1.4)这个 −1,2,−1-1,2,-1−1,2,−1 的逆矩阵被引进是因为它的特征值的形式是 2−2cosθ2-2\cos\theta2−2cosθ,所以我们可以得到 AATAA^TAAT 的 λ′s\lambda'sλ′s 和 AAA 的 σ′s\sigma'sσ′s:λ=12−2cosθ=14sin2(θ/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)λ=2−2cosθ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(2n−1)π(n=4时包括θ=93π,此时2sin2θ=1)这个特殊案例得到了当 n=4n=4n=4 时 AATAA^TAAT 的一个特征值 λ=1\lambda=1λ=1,所以 σ=λ=1\sigma=\sqrt\lambda=1σ=λ=1 时 AAA 的一个奇异值,当向量 u=(1,1,0,−1)\boldsymbol u=(1,1,0,-1)u=(1,1,0,−1) 时,有 AATu=uAA^T\boldsymbol u=\boldsymbol uAATu=u,这是一个真正的特例。
重点是画出 AAA 的 nnn 个奇异值的图形,这些数字是逐渐减小的(不像 AAA 的特征值,都是 111),但是减小的不是很陡峭。所以 SVD 只是适度的压缩这个三角形旗帜。对于希尔伯特矩阵压缩效果非常好。
四、例题
【例5】MATLAB 指令 A = rand(20, 40) 和 B = randn(20, 40) 会生成 20×40 的随机矩阵:AAA 中的元素是 000 到 111 之间均匀分布的;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;
运行结果如下:
![]() | ![]() |
---|