图片压缩：SVD（奇异值分解）之近似矩阵

写这篇博客的初衷的是想让读者浅尝一下SVD的魅力。 处于大数据时代，真的是，“SVD在手，天下我有”。最关键的是，SVD很容易理解，简单的线性代数知识就够了。

矩阵的秩

对于一个矩阵，比方如下：
$$\mathbf{A}={\left[\begin{array}{ccccc}6& 6& -2& 4& -8\\ -15& -15& 5& -10& 20\\ 12& 12& -4& 8& -16\\ 15& 15& -5& 10& -20\\ 6& 6& -2& 4& -8\end{array}\right]}_{5\times 5}$$
正常情况下，我们会以 $25$ 个字节的形式存储下来。其实呢，这个矩阵的秩为 $1$，是我随机生成的两个长度为 $5$ 向量相乘得来的：
$${\mathbf{a}}_1=[2,-5,4,5,2{]}^T,{\mathbf{a}}_2=[3,3,-1,2,-4],\mathbf{A}={\mathbf{a}}_1\times {\mathbf{a}}_2$$
所以存储起来， $10$ 个字节就够了。想象一个，如果是一个 $n\times m$ 的矩阵呢，特别是处于大数据时代的我们， $n$ 和 $m$ 可以随随便便就上万，甚至上亿，相似的方式存储这些矩阵，将会节省很多存储空间。

跟着这个思路，我们知道，对于一个大型的矩阵，如果它不是满秩，我们是可以用以上的方式来降低存储量。然而，现实生活中，这种矩阵， 比方说是每个列是某一个时刻的传感器数据，它虽然很可能不是满秩，但是也有可能接近满秩。 庆幸的是，这种大型的矩阵所包含的信息中，有一些是可以忽略的。如果我们可以求得这个矩阵的低阶近似矩阵，我们便可以达到最初的目的。

近似矩阵

为了降低存储量，我们的目标是找到一个大型矩阵的低阶近似矩阵。对于近似的理解，在向量的表示中，我们一般都是用欧几里德范数（Euclidean norm）或者欧几里德距离来衡量两个向量是否接近。在矩阵中，我们可以用弗罗贝尼乌斯范数（Frobenius norm）来衡量。所以，对于一个大型矩阵，比方说 $\mathbf{X}$，我们希望找到一个低阶的近似矩阵 $\tilde{\mathbf{X}}$。 ${\Vert \mathbf{X}-\tilde{\mathbf{X}}\Vert }_F$ 越小，两个矩阵越近似。

SVD

SVD 提供了一种系统的方法来根据主要特征确定高维数据的低阶近似。 这种技术是基于数据的，因为主要特征纯粹是从数据中发现的，不需要额外的信息。 SVD 在数值上是稳定的，并根据由数据内的主要相关性定义的新坐标系对数据进行分层。 最重要的一点是，SVD 对于任何矩阵都保证存在。

对于实数矩阵 $\mathbf{X}$ （如果是复数的话，以下的转置都是共轭转置），我们把它表示为：
$$\mathbf{X}=\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{x}}_1& {\mathbf{x}}_2& \cdots & {\mathbf{x}}_m\\ |& |& & |\end{array}\right]$$
其中，${\mathbf{x}}_i$ 是长度为 $n$ 的列向量（通常我们遇到的数据中，$n\gg m$，我们以下的表述中也基于这个假设）。存在且唯一存在一个SVD可以将矩阵 $\mathbf{X}$ 分解为：
$$\mathbf{X}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T$$
其中，$\mathbf{U}$ 和 $\mathbf{V}$ 分别是 $n\times n$ 和 $m\times m$ 的正交矩阵 （${\mathbf{U}}^T\mathbf{U}=I_n$，${\mathbf{V}}^T\mathbf{V}=I_m$）；$\mathbf{\Sigma }$ 是 $n\times m$ 的矩阵，且除了对角线上是非负数外，其余元素都是 $0$。

由于 $n⩾m$， $\mathbf{\Sigma }$ 对角线上最多有 $m$ 个非零数值，我们也可以将 $\mathbf{\Sigma }$ 表示为：
$$\mathbf{\Sigma }=\left[\begin{array}{c}\hat{\mathbf{\Sigma }}\\ 0\end{array}\right]$$
进而，我们可以将 $\mathbf{X}$ 表示为：
$$\mathbf{X}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T=[\hat{\mathbf{U}}{\hat{\mathbf{U}}}^{\perp }]\left[\begin{array}{c}\hat{\mathbf{\Sigma }}\\ 0\end{array}\right]{\mathbf{V}}^T=\hat{\mathbf{U}}\hat{\mathbf{\Sigma }}{\mathbf{V}}^T$$
也即简化的SVD（economy SVD）。值得注意的是，$\mathbf{\Sigma }$ 对角线上数值是从大到小排列的。我们将上式展开，得到：
$$\mathbf{X}=\sum _{i=1}^m{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^T={\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^T+{\sigma }_2{\mathbf{u}}_2{\mathbf{v}}_2^T+\cdots +{\sigma }_m{\mathbf{u}}_m{\mathbf{v}}_m^T$$

近似定理（Eckart-Young theorem）：
矩阵 $\mathbf{X}$ 的最优 $r$ 阶近似矩阵，${\mathbf{X}}_{\text{opt}}^r=\underset{\tilde{\mathbf{X}},\text{ s.t. rank}(\tilde{\mathbf{X}})=r}{\mathrm{argmin}}{\Vert \mathbf{X}-\tilde{\mathbf{X}}\Vert }_F$， 是截断的 $r$ 阶SVD，也即:
$${\mathbf{X}}_{\text{opt}}^r=\sum _{i=1}^r{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^T={\sigma }_1{\mathbf{u}}_1{\mathbf{v}}_1^T+{\sigma }_2{\mathbf{u}}_2{\mathbf{v}}_2^T+\cdots +{\sigma }_r{\mathbf{u}}_r{\mathbf{v}}_r^T$$

意味着，如果你期望近似矩阵为 $10$ 阶的话，那么最近似的矩阵就是：${\sum }_{i=1}^{10}{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^T$。

图片压缩

把灰度（greyscale）图片看作是一个矩阵，我想大家都很容易理解这一点。下图中，图一是一个 $1600\times 1200$ 的灰度图片（丹麦的美人鱼）。

通过上述的近似矩阵方法，我们可以分别获取 $20$ 阶、$100$ 阶 和 $250$ 阶的低阶近似矩阵（基于MATLAB）：

A = imread('./littlemermaid.jpeg');
X = double(rgb2gray(A)); %convert RGB to gray
nx = size(X,1); ny = size(X,2);

[U, S, V] = svd(X,'econ');

figure, subplot(1,4,1)
imagesc(X), axis off, colormap gray
title('Greyscale')

plotind = 2;
for r = [20 100 250] %truncation value
    Xapprox = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'; %approximate image
    subplot(1,4,plotind), plotind = plotind+1;
    imagesc(Xapprox), axis off
    title(['r=',num2str(r,'%d'),', ', num2str(100*r*(nx+ny)/(nx*ny),'%2.2f'),'% storage'])
end

对比可以看出，低阶近似矩阵需要的存储量大大降低了，而且在实际中也经常有它的用处（比方说微信可以选择发送原图或者压缩的图片，这里就可以申明使用SVD压缩。当然，也可以用其它的方法压缩）。