7.1 线性代数进行图像处理

passxgx

已于 2025-02-27 23:41:22 修改

阅读量887

点赞数 23

分类专栏： # 第7章奇异值分解（SVD）文章标签：线性代数图像处理

于 2025-02-24 22:12:55 首次发布

第7章奇异值分解（SVD）专栏收录该内容

4 篇文章

订阅专栏

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

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

二、低秩图像（例子）

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

【例1】 $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}$ $36$ 个数字变成了 $12$ 个数字，如果是 $300×300300\times300$ 个像素，那么就是 $90000$ 个数变成 $600$ 个数。如果我们提前定义好全一向量 $x\boldsymbol x$ ，我们只需要发送一个数字，这个数字就是常数灰度 $g$ ，它乘上 $xxT\boldsymbol x\boldsymbol x^T$ 得到这个矩阵。
当然第一个例子很极端，但它阐明了一个重要观点，如果能够预先定义一些特殊的向量如全一向量 $x=ones\boldsymbol x=\textrm{\pmb{ones}}$ ，那么图像处理将会变得非常高效，本质上是预选基（preselected bases）（如傅里叶基通过 FFT 进行加速）和自适应基（adaptive bases） 之间的较量。SVD 通过图像本身生成基 —— 这个是自适应的，它的计算成本可能比较高。
并不是说 SVD 总是或者通常在实际上是最有效的算法，下面这些例子的目的只是为了介绍而不是实际生产使用。

【例2】 $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}$ 这个虽说有 $3$ 种颜色，但是它的秩仍然是 $1$ ，我们仍然可以用一列乘上一行。这 $36$ 个元素甚至可能都不相同，但是只要是秩一的模式 $A=u1v1TA=\boldsymbol u_1\boldsymbol v_1^T$ 都可以压缩。但是如果秩提升至 $r = 2$ ，就需要 $u1v1T+u2v2T\boldsymbol u_1\boldsymbol v_1^T+\boldsymbol u_2\boldsymbol v_2^T$ ，例如：

【例3】内嵌的方块矩阵 $A=[1011]\boxed{A=\begin{bmatrix}1&0\\1&1\end{bmatrix}}$ 等于 $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$ 中的 $1$ 和 $0$ 可以为块状的 $1^{'} s$ 和 $0^{'} s$ ，它的秩仍为 $2$ ，我们同样只需要两项 $u1v1T+u2v2T\boldsymbol u_1\boldsymbol v_1^T+\boldsymbol u_2\boldsymbol v_2^T$ 。一个 $6×66\times6$ 的图像可以被压缩成 $24$ 个数字，一个 $N×NN\times N$ 的图像（ $N^2$ 个数字）能够压缩成来自向量 $u1,v1,u2,v2\boldsymbol u_1,\boldsymbol v_1,\boldsymbol u_2,\boldsymbol v_2$ 这四个向量的 $4 N$ 个数字。
这里 $u′s\boldsymbol u's$ 和 $v′s\boldsymbol v's$ 是最佳的选择吗？这个并不是 SVD 的选择！注意到 $u1=(1,1)\boldsymbol u_1=(1,1)$ 和 $u2=(1,0)\boldsymbol u_2=(1,0)$ 并不是正交的， $v1=(1,1)\boldsymbol v_1=(1,1)$ 和 $v2=(0,1)\boldsymbol v_2=(0,1)$ 也不正交。理论表明：正交性能够使得第二项 $c2u2v2Tc_2\boldsymbol u_2\boldsymbol v_2^T$ 更小。（SVD 按照重要性的顺序选择秩一项。）
如果 $A$ 的秩远大于 $2$ ，就比如实际图像，那么 $A$ 就会是很多秩一项相加，我们希望最小的项足够小 —— 它们可以被丢弃而不会影响视觉的质量，这时图像压缩就变为有损的了，但是好的图像压缩人的视觉几乎是无法察觉的。
问题变成了：SVD 的正交基要如何选择？

三、SVD 的特征向量

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

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)$

在例 3 中（嵌入的方块矩阵），下面是对称矩阵 $AA^T$ 和 $A^TA$ ： $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}$ 它们的特征值都是 $1$ ，所以 $λ1λ2=1\lambda_1\lambda_2=1$ ，它们的迹（对角元素之和）都是 $3$ ： $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}$ $λ1\lambda_1$ 和 $λ2\lambda_2$ 的平方根是 $σ1=5+12\sigma_1=\displaystyle\frac{\sqrt5+1}{2}$ 和 $σ2=5−12\sigma_2=\displaystyle\frac{\sqrt5-1}{2}$ 且 $σ1σ2=1\sigma_1\sigma_2=1$ 。最接近 $A$ 的秩一矩阵是 $σ1u1v1T\sigma_1\boldsymbol u_1\boldsymbol v_1^T$ ，误差仅为 $σ2≈0.6\sigma_2\approx0.6$ ，这个是理论最小值。
$AA^T$ 和 $A^TA$ 标准正交特征向量是 $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)$ 实际生活中这些计算都是由计算机完成的！（MATLAB 中可以使用 svd(A) 进行奇异值分解。）我们可验证一下 $σ1u1v1T+σ2u2v2T\sigma_1\boldsymbol u_1\boldsymbol v_1^T+\sigma_2\boldsymbol u_2\boldsymbol v_2^T$ 是否可以正确的重构矩阵 $A$ ：

$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[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)$

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

【例4】假设旗帜是两种不同颜色的三角形，左下角的三角形都是 $1$ ，右上角的三角形都是 $0$ ，主对角线都是 $1$ 。下面是当 $n = 4$ 时的图像矩阵，它是满秩的 $r = 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$ 有一组 $n$ 个奇异值 $σ\sigma$ （都是正数），SVD 会产生 $n$ 个秩一项 $σiuiviT\sigma_i\boldsymbol u_i\boldsymbol v_i^T$ ，完美的重构需要全部的 $n$ 项。
在压缩过程中小的 $σ\sigma$ 可以被丢弃，它们不会明显影响图片的质量。我们要理解并画出当 $n = 4$ 时的这些 $σ\sigma$ 以及更大的 $n$ 。注意例 3 是例 4 这个三角形矩阵当 $n = 2$ 时的特例。
下面我们手动计算，从 $AA^T$ 开始（计算机处理会不同）： $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)$ 这个 $- 1, 2, - 1$ 的逆矩阵被引进是因为它的特征值的形式是 $2−2cos⁡θ2-2\cos\theta$ ，所以我们可以得到 $AA^T$ 的 $λ′s\lambda's$ 和 $A$ 的 $σ′s\sigma'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)$ $n$ 个不同的角度 $θ\theta$ 是 $0$ ~ $π$ 之间等距离的，注意这个三对角矩阵最后一个元素是 $1$ 而不是 $2$ ，此时等距的点有 $2 n$ 个，若为 $2$ 时等距点有 $n$ 个，这个例子非常特殊： $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)$ 这个特殊案例得到了当 $n = 4$ 时 $AA^T$ 的一个特征值 $λ=1\lambda=1$ ，所以 $σ=λ=1\sigma=\sqrt\lambda=1$ 时 $A$ 的一个奇异值，当向量 $u=(1,1,0,−1)\boldsymbol u=(1,1,0,-1)$ 时，有 $AATu=uAA^T\boldsymbol u=\boldsymbol u$ ，这是一个真正的特例。
重点是画出 $A$ 的 $n$ 个奇异值的图形，这些数字是逐渐减小的（不像 $A$ 的特征值，都是 $1$ ），但是减小的不是很陡峭。所以 SVD 只是适度的压缩这个三角形旗帜。对于希尔伯特矩阵压缩效果非常好。

在这里插入图片描述

四、例题

【例5】MATLAB 指令 A = rand(20, 40) 和 B = randn(20, 40) 会生成 20×40 的随机矩阵： $A$ 中的元素是 $0$ 到 $1$ 之间均匀分布的；B 中的元素是 “钟形” 正态分布。使用 MATLAB 中的 svd 命令，画出奇异值 $σ_1$ 到 $σ_{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;

运行结果如下：