k 阶奇异值分解之图像近似

我们都知道，一般情况下，一张图像在计算机中的存储格式是三个矩阵（RGB 格式），当然也有四个矩阵（RGBA 格式）或者一个矩阵（灰度图）的情形。然而，进行数据传输的过程中如果直接从发送方把数据原封不动的传给接收方会非常浪费传输带宽，传输时延也会随之增加。在不改变通信条件的情况下，要想减少带宽占用和传输时延，只能对数据进行压缩。稍微想一下，对图像的压缩不就是对矩阵的压缩吗？矩阵压缩有很多种方法，在这里我采用 k 阶奇异值分解方法。

需求

我们需要使用 k 阶奇异值分解方法对一张图片做近似，在讲怎么对图片做近似之前，首先需要限制一些额外的条件：

1. 图片格式采用灰度图。其实采用其他格式逻辑上也是差不多的，就是矩阵的个数有所不同而已。
2. k 分别取(1,3,5,10,20,50,100)。
3. 图像处理框架我分别选择 pillow 和 scikit-image，主要是想做个对比，找到最快的方法。
4. 奇异值分解的包总共有 4 个：numpy，scipy，tensorflow 和 pytorch。其中 tensorflow 和 pytorch 既可以在 CPU 上运行，也可以在 GPU 上运行。这里考虑所有可能，做个对比，找出最快的方法！

综上所述，总共需要考虑 2✖(1+1+1*2+1*2)=12 种可能。

准备工作

准备工作总共有 2 件事，奇异值分解的原理和彩色图转灰度图（可选），我们首先看一下彩色图怎么转到灰度图。

彩色图转灰度图

如果有灰度图的话，彩色图转灰度图这一步的操作可以不进行，直接跳到奇异值分解的原理即可。这里使用两种框架：pillow 和 scikit-image，当然也有其他的图像处理框架，比如 opencv-python，只不过我没有安装过 opencv-python，所以就不去使用这个框架了，当然有这个框架的人可以自己试试。言归正传，大致逻辑：读取图片，转换图片，保存图片。我们首先看一下 pillow 框架怎么实现这个功能。

01

pillow 框架实现

pillow 框架实现这个功能非常简单，导包之后一行代码即可：

from PIL import Image
Image.open('原图.bmp').convert('L').save('黑白原图.bmp')

其中 open 方法读取图片，参数显然是图片文件的名字，返回一个图片对象，convert 方法对图像进行转换，参数对应转换方式，这里参数是 L，表示转换成灰度图，返回转换后的图片对象，save 方法用于保存图片，参数依旧是文件名。这里一行代码确实实现了图片的读取转换保存，但是，我们一会需要用 k 阶奇异值分解方法对灰度图片进行近似，在此之前，我们需要拿出灰度图的对应矩阵，如果不修改上面的代码，就需要对黑白原图.bmp 重新读取。如果这样的话，我们稍微想一下，是不是需要进行 2 次读取磁盘，一次写入磁盘，总共进行了 3 次磁盘读写。对上面的代码稍微变一下，就可以减少一次磁盘读取，我们用一个指针指向转换后的图片对象，这样，对这个指针既可以调用 save 方法，也可以通过作为 np.array 函数的参数的方式拿到图像对应的矩阵。修改后的代码如下所示。

from PIL import Image
import numpy as np
A = Image.open('原图.bmp').convert('L')
A.save('黑白原图.bmp')
A = np.array(A)

在这里我并没有通过读取磁盘中的黑白原图.bmp 的方式构造灰度图对应的图片对象，而是直接利用调用 convert 方法返回的灰度图对象，通过这种方式减少了一次磁盘的读取。在这里还有最后一件事，怎么通过一个数组创建一个图片对象？很简单，使用 Image.fromarray 方法即可，该方法有两个参数，第一个参数是数组对象，第二个参数是图片格式（和 convert 方法的参数一样）。

01

scikit-image 框架实现

scikit-image 框架实现这个功能也是非常的简单，代码如下所示：

from skimage import io, color
A = (color.rgb2gray(io.imread('原图.bmp'))*255).astype('uint8')
io.imsave('黑白原图.bmp', A)

其中 io.imread 方法用于读取图片，参数显然是图片文件的名字，需要注意的是，该方法不是返回一个图片对象，而是一个图片对象对应的一个或者多个矩阵，因此没有必要使用 np.array 函数，直接把它当成数组就行了。然后是把这个数组作为 color.rgb2gray 方法的参数调用 color.rgb2gray 方法，返回值需要注意一下，它返回的是规范化之后的灰度图矩阵，也就是说矩阵中的每个元素都是区间[0,1]的浮点数，不是 8 位二进制无符号整数。要想把矩阵的元素变成 8 位无符号整数，千万不能直接用 astype 方法转换！因为 astype 方法转换的时候类似于 C 语言的强制类型转换，会把小数点后的值抹掉。至于怎么去做，我们先反过来考虑，一个元素的值是 8 位二进制无符号整数，如何让其位于区间[0,1]内？这不就是让我手工实现 0-1 标准化吗？假设该矩阵为 A，其中某个元素值为 a，转换后的对应元素值为 b，直接利用公式：b=(a-A.min())/(A.max()-A.min())，因为 A 中的元素类型是 8 位二进制无符号整数，所以 A.min()=0，A.max()=255，代入可得 b=a/255，可是我们已知 b，需要求出 a，这太简单了，等式两边同乘 255 就可以求出 a 了，最后进行类型转换就 OK 了。在这里我利用 numpy 数组的广播机制，直接对一个数组✖255 的方式来表示对数组中每个元素✖255。然后通过调用 astype 方法进行类型转换，其参数为需要转换的数据类型。然后通过调用 io.imsave 方法进行图片的保存，该方法第一个参数是文件名，第二个参数是图像对应的数组。

奇异值分解的原理

接下来我们看到奇异值分解的原理，奇异值分解就是把一个矩阵分解成三个矩阵，公式：A=U∑V'，A 是需要进行分解的矩阵，U、∑、V 也都是矩阵，V'表示 V 的转置。

现在关键是求出 U、∑、V，我们首先计算 AA'的特征值和特征向量，然后对特征向量进行正交化和单位化，把特征值进行降序排序（特征向量也要一起带上），设排好序的非 0 特征值为λ1，λ2，……，λr，对应的特征向量为u1，u2，……，ur，令 U=(u1,u2,……,ur)，此时 U 已经出来了。然后去计算 A'A 的特征值和特征向量，然后对特征向量进行正交化和单位化，把特征值进行降序排序（特征向量也要一起带上），设排好序的非 0 特征值为μ1，μ2，……，μr，对应的特征向量为v1，v2，……，vr，令 V=(v1,v2,……,vr)，此时 V 已经出来了，如果学过线性代数的话早就可以发现 λi=μi。此时还差一个∑，∑=diag(sqrt(λ1),sqrt(λ2),……,sqrt(λr))，diag 表示主对角线是对应的元素，其余位置全为 0，sqrt 表示开算术平方根。上面是奇异值分解的数学原理，目前还差一个 k 阶该怎么解释，k 阶其实很简单，就是把特征值和特征向量提取前 k 个，当然公式也要发生一点变化：A≈U∑V'，剩下的逻辑和上面完全一样。下面我就来举个例子说一下它为什么可以压缩一个矩阵，假设有一个 1432 行 1910 列的一个矩阵，其中共有 1432✖1910=2735120 个元素，如果令 k=5 进行 k 阶奇异值分解，我们就只需要 5✖(1432+1910+1)=16715 个元素。

奇异值分解的实现

接着我们看到奇异值分解的实现，在这里我使用 6 种方法来实现：numpy、scipy、tensorflow（CPU）、tensorflow（GPU）、pytorch（CPU）、pytorch（GPU）。在此之前我们先看一下黑白原图的样子，如图所示。

01

numpy 实现

numpy 实现奇异值分解的代码很简单，如下所示：

def svd_numpy(a, k0):
    u, s, vh = np.linalg.svd(a)
    return ((u[:, :k0]*s[:k0])@vh[:k0]).astype('uint8')

在这里，我定义了一个函数，该函数有两个参数，第一个参数 a 是需要进行分解的矩阵，第二个参数 k0 是 k 阶的 k 值。通过把参数 a 作为 np.linalg.svd 的参数来调用该函数返回三个值，第一个值是矩阵 U（二维数组），第二个值是∑的对角线元素（一维数组），第三个值是 V'（二维数组），然后对矩阵做近似，返回近似之后的结果。

我们尝试令 k 分别取(1,3,5,10,20,50,100)，并采用上述两种图像处理框架，来看一下该函数的效果，并且计算一下总处理时长。numpy+pillow 的代码如下：

from time import time
from PIL import Image
import numpy as np


def svd_numpy(a, k0):
    u, s, vh = np.linalg.svd(a)
    return ((u[:, :k0]*s[:k0])@vh[:k0]).astype('uint8')


if __name__ == '__main__':
    start = time()
    A = Image.open('原图.bmp').convert('L')
    A.save('黑白原图.bmp')
    A = np.array(A)
    for k in (1, 3, 5, 10, 20, 50, 100):
        B = svd_numpy(A, k)
        Image.fromarray(B, 'L').save(f'k={k}.bmp')
    print(time()-start)

此时我们可以发现同目录下有很多 k=xxx.bmp 的图片文件，我们来看看 k=100.bmp 到底近似到什么样子了，如图所示。

和上面的原图相比还是有点模糊，下面来看一下运行时间，如图所示。

运行时间：1.2918944358825684 秒。

因为效果都是差不多的，所以就不重复给出，下面来看看 numpy+sckit-image 的代码和运行时间。

from time import time
from skimage import color, io
import numpy as np


def svd_numpy(a, k0):
    u, s, vh = np.linalg.svd(a)
    return ((u[:, :k0]*s[:k0])@vh[:k0]).astype('uint8')


if __name__ == '__main__':
    start = time()
    A = (color.rgb2gray(io.imread('原图.bmp'))*255).astype('uint8')
    io.imsave('黑白原图.bmp', A)
    for k in (1, 3, 5, 10, 20, 50, 100):
        B = svd_numpy(A, k)
        io.imsave(f'k={k}.bmp', B)
    print(time()-start)

运行时间：3.387892246246338 秒。

02

scipy 实现

scipy 实现和 numpy 几乎完全一样，只需要把上面代码的 import numpy as np 后面加上 import scipy.linalg，u, s, vh = np.linalg.svd(a)改成 u, s, vh = scipy.linalg.svd(a)就彻底地 OK 了，下面我们直接来看一下 pillow+scipy 和 scikit-image+scipy 的运行时间。

pillow+scipy 运行时间：0.3308093547821045 秒。

scikit-image+scipy 运行时间：2.4134678840637207 秒。

03

tensorflow 实现

tensorflow 实现需要考虑 CPU 和 GPU 两种情况，我们先看一下代码：

def svd_tensorflow(a, k0, device='/device:cpu:0'):
    with tf.device(device):
        a = tf.constant(a, dtype='float')
        s, u, v = tf.linalg.svd(a)
    return np.array((u[:, :k0]*s[:k0])@tf.transpose(v)[:k0], dtype='uint8')

在这里我通过 device 参数指定设备名，其他参数同上，返回值首先需要注意顺序，其次需要注意最后一个返回值返回的是矩阵 V，不是 V'。我们首先来看一下 pillow+tensorflow（CPU）的代码：

from time import time
from PIL import Image
import numpy as np
import tensorflow as tf


def svd_tensorflow(a, k0, device='/device:cpu:0'):
    with tf.device(device):
        a = tf.constant(a, dtype='float')
        s, u, v = tf.linalg.svd(a)
    return np.array((u[:, :k0]*s[:k0])@tf.transpose(v)[:k0], dtype='uint8')


if __name__ == '__main__':
    start = time()
    A = Image.open('原图.bmp').convert('L')
    A.save('黑白原图.bmp')
    A = np.array(A)
    for k in (1, 3, 5, 10, 20, 50, 100):
        B = svd_tensorflow(A, k, '/device:cpu:0')
        Image.fromarray(B, 'L').save(f'k={k}.bmp')
    print(time()-start)

我们看一下其运行时间，如图所示。

pillow+tensorflow（CPU）运行时间：1.835120677947998 秒

把运行设备换成 GPU 很简单，我们只需要把 B = svd_tensorflow(A, k, '/device:cpu:0')改成 B = svd_tensorflow(A, k, '/device:gpu:0')就 OK 了，直接看 pillow+tensorflow（GPU）的运行时间。

pillow+tensorflow（GPU）运行时间：3.783994436264038 秒。

下面我们把 pillow 改成 scikit-image，其他不变，代码的修改就参考上面的 numpy 实现的部分，照葫芦画瓢即可。我们可以首先看一下 scikit-image+tensorflow（CPU）的运行时间。

scikit-image+tensorflow（CPU）的运行时间：3.2006115913391113 秒。

下面看一下 scikit-image+tensorflow（GPU）的运行时间。

scikit-image+tensorflow（GPU）的运行时间：5.34347939491272秒。

04

pytorch 实现

pytorch 实现和 tensorflow 差不多，直接给出代码：

def svd_torch(a, k0, device='cpu'):
    if'cuda' in device:
        assert torch.cuda.is_available()
    a = torch.tensor(a, dtype=torch.float, device=device)
    u, s, v = torch.svd(a)
    a = ((u[:, :k0]*s[:k0])@v.T[:k0]).cpu()if'cuda' in device else (u[:, :k0]*s[:k0])@v.T[:k0]
    return np.array(a, dtype='uint8')

参数 device 表示设备名，不过在这里不需要像 tensorflow 那么麻烦，直接用'cpu'就可以表示使用 CPU 运行代码，但是，使用 GPU 运行代码不用'gpu'，而应该使用'cuda'。返回值的顺序和 numpy 是一样的，唯一的区别就是最后一个返回值是 V，不是 V'。

pillow+pytorch 的代码只需要根据 pillow+tensorflow 修改即可，把 import tensorflow as tf 改成 import torch，函数定义和调用以及设备名直接狸猫换太子。下面首先直接看 pillow+pytorch（CPU）运行时间。

pillow+pytorch（CPU）运行时间：0.3148210048675537 秒。

然后我们看一下 pillow+pytorch（GPU）运行时间。

pillow+pytorch（GPU）运行时间：3.394061326980591 秒。

然后我们看一下，scikit-image+pytorch（CPU）运行时间。

scikit-image+pytorch（CPU）运行时间：1.7938365936279297 秒。

最后我们看一下 scikit-image+pytorch（GPU）运行时间。

scikit-image+pytorch（GPU）运行时间：4.726118326187134 秒。

总结

从上表中可以得出以下结论：

1. pillow 的图像处理性能比 scikit-image 好。
2. 对于 tensorflow 和 pytorch 来说，使用 CPU 运行时间比使用 GPU 运行时间短，可能是因为最后转为 numpy 数组的时候需要把数据从 GPU 的显存中复制到内存中花费时间。
3. pillow+scipy 或者 pillow+pytorch（CPU）这两种组合处理速度最快，虽然在我这里 pillow+pytorch（CPU）比另一个快了一点，但是这么一点可以看成是一个误差。