21. BI - SVD 矩阵分解的实际案例：利用 SVD 进行图像压缩

本文为 「茶桁的 AI 秘籍 - BI 篇 第 21 篇」

[TOC]

Hi, 你好。我是茶桁。

上一节课的内容中，咱们学习了 SVD 矩阵分解的原理，并在最后提到了，矩阵其实是做运算的一个根基。这一节课，咱们就来举一个简单的示例，拿图片来举例。

图片本身就是一个矩阵，由长和宽来进行组成。还记得咱们在深度学习基础课上讲过卷积吗？那个时候咱们接触过图片矩阵，还有卷积核也是一个矩阵，对吧？有兴趣的可以回头去看看：(《29. 深度学习进阶 - 卷积的原理》)[https://mp.weixin.qq.com/s/9p-Fu9zgM0m1JdaPjE4Aag]。

我们现在用一个简单的方式，用一个 bmp 位图，位图只有一个通道，只需要在一个长和宽的面积上面去有一些颜色的表达。

那既然是一个矩阵了，我们就可以用奇异值把它拆成几个部分，P、S、Q 这三个部分组成。这几个分量里面的 S 会有很多的特征，这些特征里面现在假设不想要所有特征，只取前 k 个，后面那些特征把它去掉。这样我们就可以做一个算法，来对图像做压缩。我们可以看一看能不能用较小的一些存储的空间去尽可能的还原图像。

那么我们将会分为以下几个步骤：

• Step1, 将图片转为矩阵
• Step2, 对矩阵进行奇异值分解，得到 P、S、Q
• Step3, 包括特征值矩阵中的 K 个最大特征值，其余特征值设为 0
• Step4, 通过 P, S', Q 得到新的矩阵 A'， 对比 A'与 A 的差别。

我平时很喜欢用 iPhone 拍照，平时拍照的时候大家有没有关注一下自己拍摄的照片大概是多大？一般来说，基本是在 2-4M 左右吧？其实我们传输的时候，特别是用微信传输，大家应该能看到一个「发送原图」的选项，也就是说，我们发送的那个图片如果没有勾选这个选项，都是经过压缩的。

你把这个图像另保存在电脑里，会发现它不是原图格式，但是也可以看的比较清楚，虽然信息有一定的损失，大小大约是在 200K 左右。从 4 兆到 200K，我们的信息其实只有百分之五左右，下降了大概 20 倍。也就说他用了 5%的信息，但是能还原出来绝大部分的一些内容。那你有没有好奇这个技术是怎么做到的？而且这个技术在图像里是个通用技术。

今天我就带着大家用 SVD 做一版图像压缩的算法，看一看它能不能用很小的一些信息帮我们保存尽可能多的一些内容。

SVD 它有一个价值就是把特征抽取出来，而且还把特征的权重从大到小做了个排序。我们可以通过中间 S 这个矩阵，就可以在对角线上看到它。

我们现在的想法就是把一些不太想要的后面那些特征给它设置为 0，只提一些关键特征。有了关键特征，我们可以把这些图像再去做一些还原，还原出来的这些图像跟原图之间做一个对比。

这次用的是一张街拍的生活照，是我从同事朋友圈 down 的，希望她不会怪我。在拿到图像后我做了一些处理，将其处理为灰度，因为这次我们需要用到一个单通道的图像。

大家应该都理解通道的概念吧？JPG 是我们最常见的格式，我们可以看到四个通道，R、G、B 以及合成通道，其实严格意义上来说，一张 JPG 只包含三个标准通道，R、G、B，三个通道任意关闭一个，合成通道都是无效的。那 R、G、B 表示的就是红，绿，蓝三个颜色。

现在为了方便起见，我们不去做三个通道的图像，将图像转为灰度之后输出为 BMP，8 位图。

我们来看原始图像

当然，这是我用 plt show 出来的，并且大家还是自己去找图片去做测试，这张图就恕我不提供了。这是图片信息：

我们需要将图片先读取进来，然后使用 NumPy 将其转为矩阵赋值给 A。

# 加载图片
image = Image.open('assets/1221-1202.bmp')
A = np.array(image)

然后我还做了个展示，看看提取出来的 A 是否可以正常显示，显示结果就如我上面贴图一样。

# 显示原图
plt.imshow(A, cmap=plt.cm.gray, interpolation='nearest')
plt.show()

有了原始图像以后，然后我们该做什么了？当然是拆矩阵对吧？要将这个图像矩阵拆成三块, P, lambda 和 Q：

# 对图像矩阵 A 进行奇异值分解，得到 P、S、Q
p, s, q = svd(A, full_matrices=False)

接着，我们现在做法就是抽它的特征。这里，咱们写一个函数，用于从 S 里面抽取几个关键特征。

# 取前 K 个特征，对图像进行还原
def get_image_feature(s, k):
  # 对于 S，值保留前 k 个特征值
  s_temp = np.zeros(s.shape[0])
  s_temp[0:k] = s[0:k]
  s = s_temp * np.identity(s.shape[0])

接着，需要对这个函数进行补全，我们不仅希望它提取特征，对于提取后的特征还原一个 temp 新矩阵，然后将它显示出来。

def get_image_feature(s, k):
  ...

  # 用新的 s_temp, 以及 p，q 重构 A
  temp = np.dot(p, s)
  temp = np.dot(temp, q)
  plt.imshow(temp, cmap=plt.cm.gray, interpolation='nearest')
  plt.show()

接着我们来传参调用函数，参数包括 s 和 k，s 是代表了特征从大到小的一个顺序关系。我用前 5 个来进行提取和还原：

get_image_feature(s, 5)

然后 50 个

get_image_feature(s, 50)

在接着是 500 个

get_image_feature(s, 500)

下面来给大家分析一下这个函数，首先，我们得到一个全 0 的 zeros

s_temp = np.zeros(s.shape[0])

它把它所有的零都设置上，打印出来应该是这样一个矩阵：

[0. 0. 0. ... 0. 0. 0.]

然后又把这个 s 里的 0 到 k 给它还原出来, k 是我们传参传进来的。

s_temp[0:k] = s[0:k]

这样我们可以想象一下，只有前面这几个值是有价值的，后面都为 0 了。

然后得到的这个有价值的矩阵 s_temp 乘上一个 identity，identity 就是我们上节课讲到的单位矩阵，称之为 i，对角线为 1 的矩阵叫做单位矩阵。

s = s_temp * np.identity(s.shape[0])

我们有了 s、p 和 q，p 和 q 是原来拆出来的内容，是不会发生变化的。这个 s 乘上 p 得到一个值赋值给一个临时值 temp，再拿得到的结果和 q 进行相乘，继续重新赋值 temp

temp = np.dot(p, s)
temp = np.dot(temp, q)

这样我们就会近似还原一张图，接着要做的事情就是跟我们展示最开始的 A 矩阵一样，将 temp 给展示出来就可以了：

plt.imshow(temp, cmap=plt.cm.gray, interpolation='nearest')
plt.show()

我们回头去看看最后展示出来的那张图，就是 k=500 的时候的那张图，前两个明显有差距我们不需要仔细看。k=500 那张图仔细看，和原图还是有区别的。最明显的，阴影部分的深度没有那么大，对吧？那有可能是在输出色阶上范围低了一点，这部分特征是丢失了。

那么，矩阵对角线里面是有权重的，它的特征值个数一般是多少？这个要看向量的维度。我之前展示了这张图片的信息，我们知道这张图片现在向量应该是 3840*2160，所以它的特征维度我们来猜一猜一般会是多少维？

这里全部的特征个数应该有一个规律，实际上是应该小于等于(3840, 2160)的最小值，也就是小于等于 min(3840, 2160)。也就是你长和宽里面的最小值要比它小，有可能就是 2160，所以它的上限是 2160，最多有可能是 2000 多维。通常情况下很有可能就是 2160。

现在 2,000 多维里面我们只取了 5 个维度, 我们脑海里过滤一下，5 个维度会不会这个图像就花掉了？图像原本是由 2,000 多个维度合并而成的，现在我们只要从大到小的前五个，事实也是如此。之前展示的 k=5 的图像确实是花的，能看出来原图是什么样子吗？前五个的信息量其实已经还挺大的，但是你用的数据太少，基本上是不可能的。你想 5 除上 2,000，这只有多少的压缩空间啊？1%都不到。

再看一看 50 个维度，原来是 2160 个维度，现在变成了 50，信息其实也是用的非常少。但是这 50 个特征，已经能看出来原图的样子，有一点眉目了。这 50 个特征应该基本上就能看出来原来图像的一个概念。所以我们用这个方式就可以很好的帮你来做还原。

50 个可以看到了，如果我们用更多的 500 个，可以看到由原来很模糊到现在很清楚，这个差距其实是非常明显的。那以上就是 SVD 的进行图片压缩的过程。

所以 SVD 可以帮你来做一个降维的处理。这里的降维我们先说结论，可以用很少的信息，大概 10%左右就可以还原大部分的一些信息内容，信息可以还原出来 90%左右。

那为什么是这样的一个比值呢？我们以刚才的 k 等于 50 为例，k 等于 50 是怎么保存的。

(m + 1 + n) * k = ?

我们的 m 乘上 k 是前面一个矩阵，对于后面那些都为 0 的那些部分我们的存储空间是不需要存储的，因为它是乘法，是没有意义的。这是 m 乘上 k，中间这个部分的应该是单位的对角阵，单位对角阵的话现在应该大小应该是 1 乘上 k。因为只要把它保存向量就好了，就像刚才我们看到那个 s 是一样的。后面这个部分的应该就是 n，那就是 k 乘上 n，就是 2160 乘上 50 这个信息，也就是 Q 里面只需要存2160*50。

前面只需要存3840*50，中间只需要存 50 个，后面是2160*50，我把所有的元素的个数都给它存出来，这些是你一定要存的信息，(3840 + 1 + 2160) * 50 = 300050，大概有 300050 个元素。看起来虽然很多，但是相比原来这个信息量只有多少呢？我们来看看原图信息量：3840 * 2160 = 8294400，我们做一下对比，300050/(3840 * 2160) = 0.036175...， 大概是不到 4%的一个占比。

所以说我们其实可以只用了差不多 10%的信息可以还原出来差不多 90%以上的信息。可以对比一下，以上就把一个图像的压缩的原理简单的给大家讲完了。背后使用的工具是 SVD，SVD 可以很好的帮我们来分析一个矩阵中哪些成分是关键的，哪些成分不是关键的。这样我们就可以对一个矩阵去提取它的关键特征来做一些还原。

这个例子中，我们是先用到了一个图像压缩领域中让你去了解如何提取图像中的关键特征，又如何把一些重要的 top-k 特征做了一个近似还原。结论就是，10%的信息可以相当于 90%的信息量。

我们来看，这个技术是不是感觉起来挺神奇的？那想象一下，我们在微信里面传的原始图像，原来 4 兆多，虽然只用了 200K，但同样可以得到很清晰的一个图像，就是以上的一个原理。

大家在课后，可以拉取我的代码去跑一跑，不过需要换成你们自己的图片了，这张图片就不提供了。大家可以自己去体验一下，写一个图像的压缩工具。