7.1 线性代数进行图像处理

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

$A$ 的奇异值定理就是 $A^{T}A$ 和 $A{A}^T$ 的特征值定理。 这是本章的内容的预览。$A$ 有两个奇异向量的集合（$A^{T}A$ 和 $A{A}^T$ 的特征向量），有一个是正奇异值的集合（因为 $A^{T}A$ 和 $A{A}^T$ 有相同的正特征值）。通常 $A$ 是长方形的矩阵，但是 $A^{T}A$ 和 $A{A}^T$ 是方阵，对称的且至少是半正定的。
奇异值分解（The Singular Value Decomposition：SVD）可以将任意的矩阵分解成简单的片（pieces）。 每个片都是一个列向量乘一个行向量，一个 $m\times n$ 的矩阵有 $m\times n$ 个元素（如果表示的是一个图像则会是很大的数字），但是列和行仅有 $m+n$ 个分量，远少于 $m\times n$。这些列乘行 (colum)(row) 的小片是可以非常快速处理的全尺寸大小矩阵（full size matrices）—— 它们仅需要 $m+n$ 个数。
不寻常的是，SVD 的图像处理应用是在它所依赖的代数矩阵之前出现的，下面会从只包含一个或两个片的简单图像开始。现在我们将图像想象成一个大的长方形矩阵，元素 $a_{ij}$ 表示的是图像所有像素的灰度值，可以将像素想象成一个小的正方形，从左下角开始向右 $i$ 步和向上 $j$ 步，它的灰度值是一个数字（通常的范围是 $0\le a_{ij}\le 256=2^8$），一个纯白的像素是 $a_{ij}=255=11111111$，计算机将 $255$ 用二进制表示时，这个数字就是 $8$ 个 $1$。
目前我们已经知道一个图像有 $m\times n$ 个像素，每个像素使用 $8$ 位（$0$ 或 $1$）表示它的灰度，变成了一个 $m\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=\left[\begin{array}{cccccc}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{array}\right]而是发送A=\left[\begin{array}{c}1\\ 1\\ 1\\ 1\\ 1\\ 1\end{array}\right]\left[\begin{array}{cccccc}1& 1& 1& 1& 1& 1\end{array}\right]$$$36$ 个数字变成了 $12$ 个数字，如果是 $300\times 300$ 个像素，那么就是 $90000$ 个数变成 $600$ 个数。如果我们提前定义好全一向量 $\mathbf{x}$，我们只需要发送一个数字，这个数字就是常数灰度 $g$，它乘上 $\mathbf{x}{\mathbf{x}}^T$ 得到这个矩阵。
当然第一个例子很极端，但它阐明了一个重要观点，如果能够预先定义一些特殊的向量如全一向量 $\mathbf{x}=\text{ones}$，那么图像处理将会变得非常高效，本质上是预选基（preselected bases）（如傅里叶基通过 FFT 进行加速）和自适应基（adaptive bases） 之间的较量。SVD 通过图像本身生成基 —— 这个是自适应的，它的计算成本可能比较高。
并不是说 SVD 总是或者通常在实际上是最有效的算法，下面这些例子的目的只是为了介绍而不是实际生产使用。

【例2】 $$\begin{array}{l}“\text{ace}三色旗”\\ 法国国旗A\\ 意大利国旗A\\ 德国国旗A^T\end{array}不要发送A=\left[\begin{array}{cccccc}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{array}\right]发送A=\left[\begin{array}{c}1\\ 1\\ 1\\ 1\\ 1\\ 1\end{array}\right]\left[\begin{array}{cccccc}a& a& c& c& e& e\end{array}\right]$$这个虽说有 $3$ 种颜色，但是它的秩仍然是 $1$，我们仍然可以用一列乘上一行。这 $36$ 个元素甚至可能都不相同，但是只要是秩一的模式 $A={\mathbf{u}}_1{\mathbf{v}}_1^T$ 都可以压缩。但是如果秩提升至 $r=2$，就需要 ${\mathbf{u}}_1{\mathbf{v}}_1^T+{\mathbf{u}}_2{\mathbf{v}}_2^T$，例如：

【例3】内嵌的方块矩阵 $\boxed{A=\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]}$ 等于 $\boxed{A=\left[\begin{array}{c}1\\ 1\end{array}\right]\left[\begin{array}{cc}1& 1\end{array}\right]-\left[\begin{array}{c}1\\ 0\end{array}\right]\left[\begin{array}{cc}0& 1\end{array}\right]}$

这里矩阵 $A$ 中的 $1$ 和 $0$ 可以为块状的 $1^{\prime }s$ 和 $0^{\prime }s$，它的秩仍为 $2$，我们同样只需要两项 ${\mathbf{u}}_1{\mathbf{v}}_1^T+{\mathbf{u}}_2{\mathbf{v}}_2^T$。一个 $6\times 6$ 的图像可以被压缩成 $24$ 个数字，一个 $N\times N$ 的图像（$N^2$ 个数字）能够压缩成来自向量 ${\mathbf{u}}_1,{\mathbf{v}}_1,{\mathbf{u}}_2,{\mathbf{v}}_2$ 这四个向量的 $4N$ 个数字。
这里 ${\mathbf{u}}^{\prime }s$ 和 ${\mathbf{v}}^{\prime }s$ 是最佳的选择吗？这个并不是 SVD 的选择！注意到 ${\mathbf{u}}_1=(1,1)$ 和 ${\mathbf{u}}_2=(1,0)$ 并不是正交的，${\mathbf{v}}_1=(1,1)$ 和 ${\mathbf{v}}_2=(0,1)$ 也不正交。理论表明：正交性能够使得第二项 $c_2{\mathbf{u}}_2{\mathbf{v}}_2^T$ 更小。（SVD 按照重要性的顺序选择秩一项。）
如果 $A$ 的秩远大于 $2$，就比如实际图像，那么 $A$ 就会是很多秩一项相加，我们希望最小的项足够小 —— 它们可以被丢弃而不会影响视觉的质量，这时图像压缩就变为有损的了，但是好的图像压缩人的视觉几乎是无法察觉的。
问题变成了：SVD 的正交基要如何选择？

三、SVD 的特征向量

大部分图像的特征向量并不正交，而且特征向量 ${\mathbf{x}}_1,{\mathbf{x}}_2$ 只提供一组向量，我们想要两组向量（${\mathbf{u}}^{\prime }s$ 和 ${\mathbf{v}}^{\prime }s$）。这两个问题的解决方法正是 SVD 的思想：
使用 $A{A}^T$ 的特征向量 $\mathbf{u}$ 和 $A^{T}A$ 的特征向量 $\mathbf{v}$。
因为 $A{A}^T$ 和 $A^{T}A$ 自动就是对称的（通常并不相等！），所以特征向量 ${\mathbf{u}}^{\prime }s$ 是一组正交向量，特征向量 ${\mathbf{v}}^{\prime }s$ 是另一组正交向量，我们可以并且需要将它们单位化：$||{\mathbf{u}}_i||=1$ 和 $||{\mathbf{v}}_i||=1$，这样我们的秩 $2$ 矩阵就是 $A={\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^T+{\sigma }_2{\mathbf{u}}_2{\mathbf{v}}_2^T$。其中 ${\sigma }_1$ 和 ${\sigma }_2$ 这些数字的大小将决定它们在压缩过程中是否可以被忽略。我们保留大的 ${\sigma }^{\prime }s$，丢弃小的 ${\sigma }^{\prime }s$。
SVD 中的 ${\mathbf{u}}^{\prime }s$ 称为左奇异值向量（left singular vectors），它们是 $A{A}^T$ 的单位特征向量；${\mathbf{v}}^{\prime }s$ 称为右奇异值向量（right singular vectors），它们是 $A^{T}A$ 的单位特征向量。其中 ${\sigma }^{\prime }s$ 就是奇异值（singular values），等于 $A{A}^T$ 和 $A^{T}A$ 特征值的平方根：

SVD 的选择 $A{A}^T{\mathbf{u}}_i={\sigma }_i^2{\mathbf{u}}_i{A}^{T}A{\mathbf{v}}_i={\sigma }_i^2{\mathbf{v}}_{i}A{\mathbf{v}}_i={\sigma }_i{\mathbf{u}}_i(\mathrm{7.1.1})$

在例 3 中（嵌入的方块矩阵），下面是对称矩阵 $A{A}^T$ 和 $A^{T}A$：$$A{A}^T=\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\left[\begin{array}{cc}1& 1\\ 0& 1\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 2\end{array}\right]A^{T}A=\left[\begin{array}{cc}1& 1\\ 0& 1\end{array}\right]\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]=\left[\begin{array}{cc}2& 1\\ 1& 1\end{array}\right]$$它们的特征值都是 $1$，所以 ${\lambda }_1{\lambda }_2=1$，它们的迹（对角元素之和）都是 $3$：$$\det \left[\begin{array}{cc}1-\lambda & 1\\ 1& 2-\lambda \end{array}\right]={\lambda }^2-3\lambda +1=0得到{\lambda }_1=\frac{3+\sqrt{5}}{2}和{\lambda }_2=\frac{3-\sqrt{5}}{2}$$${\lambda }_1$ 和 ${\lambda }_2$ 的平方根是 ${\sigma }_1=\frac{\sqrt{5}+1}{2}$ 和 ${\sigma }_2=\frac{\sqrt{5}-1}{2}$ 且 ${\sigma }_1{\sigma }_2=1$。最接近 $A$ 的秩一矩阵是 ${\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^T$，误差仅为 ${\sigma }_2\approx 0.6$，这个是理论最小值。
$A{A}^T$ 和 $A^{T}A$ 标准正交特征向量是$${\mathbf{u}}_1=\left[\begin{array}{c}1\\ {\sigma }_1\end{array}\right]{\mathbf{u}}_2=\left[\begin{array}{c}{\sigma }_1\\ -1\end{array}\right]{\mathbf{v}}_1=\left[\begin{array}{c}{\sigma }_1\\ 1\end{array}\right]{\mathbf{v}}_2=\left[\begin{array}{c}1\\ -{\sigma }_1\end{array}\right]全部除以\sqrt{1+{\sigma }_1^2}(\mathrm{7.1.2})$$实际生活中这些计算都是由计算机完成的！（MATLAB 中可以使用 svd(A) 进行奇异值分解。）我们可验证一下 ${\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^T+{\sigma }_2{\mathbf{u}}_2{\mathbf{v}}_2^T$ 是否可以正确的重构矩阵 $A$：

$\boxed{A=\left[\begin{array}{cc}{\mathbf{u}}_1& {\mathbf{u}}_2\end{array}\right]\left[\begin{array}{cc}{\sigma }_1& \\ & {\sigma }_2\end{array}\right]\left[\begin{array}{c}{\mathbf{v}}_1^T\\ {\mathbf{v}}_2^T\end{array}\right]}$ 或更简单的 $\boxed{A\left[\begin{array}{cc}{\mathbf{v}}_1& {\mathbf{v}}_2\end{array}\right]=\left[\begin{array}{cc}{\sigma }_1{\mathbf{u}}_1& {\sigma }_2{\mathbf{u}}_2\end{array}\right])}(\mathrm{7.1.3})$

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

【例4】假设旗帜是两种不同颜色的三角形，左下角的三角形都是 $1$，右上角的三角形都是 $0$，主对角线都是 $1$。下面是当 $n=4$ 时的图像矩阵，它是满秩的 $r=4$，所以它可逆：$$\begin{array}{l}三角形\\ 旗帜矩阵\end{array}A=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 1& 0& 0\\ 1& 1& 1& 0\\ 1& 1& 1& 1\end{array}\right]和A^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ -1& 1& 0& 0\\ 0& -1& 1& 0\\ 0& 0& -1& 1\end{array}\right]$$由于矩阵满秩，所以 $A$ 有一组 $n$ 个奇异值 $\sigma$（都是正数），SVD 会产生 $n$ 个秩一项 ${\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^T$，完美的重构需要全部的 $n$ 项。
在压缩过程中小的 $\sigma$ 可以被丢弃，它们不会明显影响图片的质量。我们要理解并画出当 $n=4$ 时的这些 $\sigma$ 以及更大的 $n$。注意例 3 是例 4 这个三角形矩阵当 $n=2$ 时的特例。
下面我们手动计算，从 $A{A}^T$ 开始（计算机处理会不同）：$$A{A}^T=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 2& 2& 2\\ 1& 2& 3& 3\\ 1& 2& 3& 4\end{array}\right]和(A{A}^T{)}^{-1}=(A^{-1}{)}^T{A}^{-1}=\left[\begin{array}{cccc}2& -1& 0& 0\\ -1& 2& -1& 0\\ 0& -1& 2& -1\\ 0& 0& -1& 1\end{array}\right](\mathrm{7.1.4})$$这个 $-1,2,-1$ 的逆矩阵被引进是因为它的特征值的形式是 $2-2\cos \theta$，所以我们可以得到 $A{A}^T$ 的 ${\lambda }^{\prime }s$ 和 $A$ 的 ${\sigma }^{\prime }s$：$$\lambda =\frac{1}{2-2\cos \theta }=\frac{1}{4{\sin }^2(\theta /2)}得到\sigma =\sqrt{\lambda }=\frac{1}{2\sin (\theta /2)}(\mathrm{7.1.5})$$$n$ 个不同的角度 $\theta$ 是 $0$ ~ $\pi$ 之间等距离的，注意这个三对角矩阵最后一个元素是 $1$ 而不是 $2$，此时等距的点有 $2n$ 个，若为 $2$ 时等距点有 $n$ 个，这个例子非常特殊：$$\theta =\frac{\pi }{2n+1},\frac{3\pi }{2n+1},\cdots ,\frac{(2n-1)\pi }{2n+1}(n=4时包括\theta =\frac{3\pi }{9}，此时2\sin \frac{\theta }{2}=1)$$这个特殊案例得到了当 $n=4$ 时 $A{A}^T$ 的一个特征值 $\lambda =1$，所以 $\sigma =\sqrt{\lambda }=1$ 时 $A$ 的一个奇异值，当向量 $\mathbf{u}=(1,1,0,-1)$ 时，有 $A{A}^T\mathbf{u}=\mathbf{u}$，这是一个真正的特例。
重点是画出 $A$ 的 $n$ 个奇异值的图形，这些数字是逐渐减小的（不像 $A$ 的特征值，都是 $1$），但是减小的不是很陡峭。所以 SVD 只是适度的压缩这个三角形旗帜。对于希尔伯特矩阵压缩效果非常好。

四、例题

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

运行结果如下：