错误使用 svd svd 的输入不能包含 nan 或 inf。_使用SRC算法实现基于投票机制的遮挡人脸图像分类实验...

使用SRC算法实现基于投票机制的遮挡人脸图像分类实验

1. 实验目的（问题描述）

该实验要求使用SRC算法实现基于投票机制的人脸遮挡分类实验。简单来说，通过无遮挡人脸训练集，来为有遮挡的人脸进行分类。SRC在此更像是一个分类的框架（方法论），而OMP(Orthogonal Matching Pursuit)是对于给定字典和样本，求样本在该字典下稀疏表达的算法，是SRC框架下的核心算法。在本次实验中，我实现了整个流程，并且做了正确性证明，给出一些问题的分析及数学原理的解释，以及尝试性的速度优化。

1. 算法描述和分析

在进行具体的编程实践之前，我首先阅读了张老师提供的[1]。它使用稀疏表达方法进行人脸识别的分类。该方法全称为Sparse Representation-based Classification ，简称SRC。该论文通过实验验证了-Minimization比-Minimization更能够确保解的稀疏性，同时也提出了稀疏表达可以用来处理分类、降噪问题。在进行程序实现之前，我使用类似自然语言来重新描述算法大框架，并且将分类问题具体化为人脸识别问题。

首先不考虑分块。

1. 使用训练样本构建矩阵，其中的也是矩阵，它是第i个人对应的矩阵，包含了14个列向量，每个列向量对应一个预处理过的训练样本。在本题中，共有1400个训练样本，所以的列数为1400。
2. 对于一个训练样本，使用OMP算法求的稀疏解。
3. 遍历所有的人，对于每一个人，经过“投影函数”（即题目中的类似“one-hot”编码，可以理解为mask），只保留在该人训练集中对应位置的值，其余值全部置零。比如，第二个人在矩阵中的列为[14-27]，那么经过投影函数后，稀疏解的1400列中，只有这14列的值被保留。
4. 使用经过投影函数投影的稀疏向量与相乘进行重构，得到（如上所述，一个人对应一个投影函数，与相乘产生一个）。计算与的残差（欧氏距离），残差最小的类即为分类结果。

在考虑分块后，对于每个块都走一遍上述的流程。也就是说，每个位置的块对应一个字典，然后通过投票来得到最终结果。

该项目并没有涉及到字典学习，而是直接把降采样的patch拉伸成纵列变成原子。因此，该项目的核心算法就是OMP算法了。为了理解OMP算法并实现，我花了不少的时间复习欠缺的线代知识，同时询问助教老师，解决了一些疑惑（感谢）！

OMP算法的数学形式如下（取自课程PPT lecture3 第15页）：

OMP算法是一种贪婪算法，它每次选择残差投影最大的向量加入（字典列向量的子集），然后重新计算表达，后重新计算残差，依次迭代。OMP先对于BP算法有很多优点，比如一定可以在特定迭代次数内收敛等。

开始之前的疑惑及解答

1. OMP算法的迭代要几轮？

答：

从第一个角度考虑，迭代次数应该等于秩的最大可能值，也就是，其中为行数，为列数。这是因为，行秩恒等于列秩，设其为，即最多可以用个列向量表达在字典的列空间中任意维的向量。所以不需要迭代1400次。迭代sample_length次即可。每次迭代产生一个非零权重，这些权重随后对应放至1400维的对应位置即可。

从第二个角度考虑，OMP中计算的本质解一个最小二乘解，即

为一个最小二乘问题。而最小二乘问题之所以是最小二乘，不能直接求得精确解，就是因为该方程组是超定的，也就是说，Ak的行数必然大于列数。因为每次迭代增加一列，所以显然迭代次数不能超过行数。

1. 为什么在OMP过程中会报奇异矩阵（Singular Matrix）不可求逆的错误？以及应对措施？

首先错误截图如下：

也就是说，为奇异矩阵，也就是不满秩矩阵，不能求逆。为什么会不满秩呢？

由线性代数知识可得，[2]，而Ak的来源A中的原子。不满秩的原因是共线性。当Ak列数比较多的时候，产生共线性实在是太正常了。比如说两个patch都是围巾，几乎都是一团黑色，就呈严重的共线性。因此必须从数学上解决这个问题，而不能希望数据比较好。经过查阅资料，伪逆（Moore-Pseudo inverse）可以较好地解决这个问题。它通过SVD进行求伪逆，虽然速度比求逆要慢一些，但是可以完美解决这个问题。

1. 归一化具体是怎么做的？为什么？

归一化的目标是使得向量的范数为，推导如下：

设向量，其范数。现在求变换,使得.设，则

也就是说，只需要把每个元素都除以范数，其范数自然会等于1.

1. 算法实现

在数据预处理过程中，需要将图片resize后分成patch，然后展平放入训练集和测试集。训练集就是字典，我使用A来标记，其shape为(9,210,1400)，是一个numpy的float数组。第一维对应着patch_index，这样可以避免创建9个变量，也方便更改。同理，测试集为TEST，其shape为(9,210,1200)。

带注释的OMP算法实现如下：

import numpy as np


def f(y, A, sample_length):
    # y=Ax
    # A: 过完备字典
    cols = []
    # K：循环次数，据分析前面分析，等于样本维数
    K = sample_length
    # Ak:行数>=列数的字典，由已经选取的原子向量构成。
    Ak = np.zeros([sample_length, 0])
    # 残差
    r = y[:, np.newaxis]
    for i in range(K):
        # 残差对每个原子向量投影，可以通过点乘字典的转置并取绝对值实现
        projections = np.abs(np.dot(A.T, r))
        # 找打投影长度最长的原子向量
        best_atom_index = np.argmax(projections)
        # 将该列索引号其加入索引号子集
        cols.append(best_atom_index)
        # 并将该列加入“子集字典”
        Ak = np.c_[Ak, A[:, best_atom_index]]
        # 重新计算表达
        xk = np.linalg.pinv(Ak.T.dot(Ak)).dot(Ak.T).dot(y[:, np.newaxis])
        # 重新计算残差
        r = y[:, np.newaxis] - Ak.dot(xk)

    # 迭代结束后，cols已经有了210个列号，xk也有了210维，cols决定着如何将xk的元素放入X的对应位置
    # X是1400维的系数向量
    X = np.zeros(shape=(1400,), dtype=np.float64)
    # 按照位置的对应关系放入非零元素
    X[cols] = xk.flatten()

return X

分块SRC算法：

我总共测试了两种降采样和分块的模式。仅附上图像大小为45*42，分块形式如下的代码。（另一种降采样和分块的模式仅仅是参数的区别）

import time
from PIL import Image
import numpy as np
import os
import cv2
from sklearn.preprocessing import normalize
from OMP import f
from collections import Counter

img_folders = r"H:JuniorYearISTUDY表达学习大项目AR"
# 读取数据集文件列表
files = os.listdir(img_folders)
# 降采样至该分辨率
down_sample_to = (42, 45)

# 计算样本长度,即字典行数
sample_length = (down_sample_to[0] // 3) * (down_sample_to[1] // 3)
print(sample_length)
# 搭建字典骨架
A = np.zeros(shape=(9, sample_length, 1400), dtype=np.float64)
# 搭建测试集骨架
TEST = np.zeros(shape=(9, sample_length, 1200), dtype=np.float64)  # 相同预处理的测试样本

# 创建字典,键是人的编号.值是该人对应的mask，该mask在字典对应atom位置为1，其余位置为0.一个mask由14个1,1386个0
mask_train = {}
# 初始化mask
for file_name in files:
    mask_train[file_name[0:5]] = []
# 构造行列分割点
row_split = [0]
col_split = [0]
for i in range(1, 4):
    row_split.append(down_sample_to[1] // 3 * i)
for i in range(1, 4):
    col_split.append(down_sample_to[0] // 3 * i)

# 指向A和TEST列标的游标,在构建A和TEST时使用
col_train = 0
col_test = 0

# 文件名中的index与训练/测试属性的对应关系
train_index = list(range(1, 8)) + list(range(14, 21))
test_index = list(range(8, 14)) + list(range(21, 27))

for file_name in files:
    img_path = os.path.join(img_folders, file_name)
    # 提取文件名中的编号
    index = int(file_name[:-4].split('-')[2])
    # 打开图像文件
    img = Image.open(img_path)
    img = np.array(img)
    # 降采样
    down_sampled_img = cv2.resize(img, down_sample_to)
    # 分块,i为某块所在的行,j为列
    for i, row_inx in enumerate(range(0, 3)):
        for j, col_inx in enumerate(range(0, 3)):
            # patch_index为块的索引,和上图对应
            patch_index = i * 3 + j
            # 根据行列分割点,分割出patch
            patch = down_sampled_img[row_split[i]:row_split[i + 1], col_split[j]:col_split[j + 1]]
            # 转换格式
            patch = patch.astype(np.float64)
            # 展平
            patch = patch.reshape(-1, 1)
            # l2归一化,即向量的每个元素都除以l2范数
            patch = normalize(patch, 'l2', axis=0)
            # 根据文件名中的index,对应放入字典或测试集
            if 1 <= index <= 7 or 14 <= index <= 20:
                A[patch_index, :, col_train] = patch.flatten()
            else:
                TEST[patch_index, :, col_test] = patch.flatten()

    # 字典或测试集的游标+1
    if index in train_index:
        mask_train[file_name[0:5]].append(col_train)
        col_train += 1
    else:
        col_test += 1

# 将所有人在字典中对应列的index转为mask（类似于one-hot）
for key in mask_train.keys():
    mask_array = np.zeros(1400, dtype=np.float64)
    mask_array[mask_train[key]] = 1
    mask_train[key] = mask_array

# 同样,指向TEST的游标
col_test = 0

# 预测正确的测试样本数目
right = 0

for person_truth in mask_train.keys():
    # 第一层for选取一个人
    # 第二层for的一次迭代选取该人对应的一个测试样本
    for i in test_index:
        start = time.time()
        # votes用于9个patch的投票
        votes = []
        # 遍历九个patch
        for patch_index in range(9):
            # 残差.键是人的编号。值是某测试样本经投影后重新表达得到的residual
            residuals = {}
            # y某张测试图片的某个patch
            y = A[patch_index, :, col_test]
            # X是OMP算法计算出来的稀疏解
            X = f(y, A[patch_index], sample_length)
            # 将X向字典里每个人的位置依次投影，每投影一次就计算一个重建后的残差。
            for guessed_person in mask_train.keys():
                projected_X = X * mask_train[guessed_person]
                r = np.linalg.norm(y[:, np.newaxis] - A[patch_index].dot(projected_X))
                residuals[guessed_person] = r
            # 残差最小的类为预测的类
            predicted = min(residuals, key=residuals.get)
            # 该类多了一票
            votes.append(predicted)
        # 列指针后移
        col_test += 1

        # 进行计数
        votes = Counter(votes)
        # person_truth为真实标签,votes.most_common(1)返回票数最高的标签.若票数最高的不止一个,则返回下标靠前的.
        print(person_truth, votes.most_common(1))
        # 如果预测正确，则预测正确数+=1
        if person_truth == votes.most_common(1)[0][0]:
            right += 1
        end = time.time()
        print(end - start)
    # 打印每个人的测试样本的准确率（12张图的准确率）
    print("准确率: ", right / 12)

1. 实验数据及数据预处理

我把文件名解析为人的编号和每个人的index。比如m-001-01.pgm，m-001就是人的编号，01就是m-001的index。

此外，还需要实现前文提到的“投影函数”。我在此使用python的字典来构造，该字典的键为人的编号，值为14为mask，如：

使用与文件名对应的编号，而不是一个递增的自然数列，降低了出错的概率。

1. 实验结果及实验分析

在评测实验结果时，使用了所有的1200个测试样本。

5.1 未分块时的测试

5.1.1 测试结果

一开始为了确保对大框架理解正确，我是没有进行分块的。我将图片降采样到12*10整体进行测试，得到正确率只有16.7%。我认为原因有以下两个：

5.1.2 原因分析

1. 主要原因：我认为和数据集有关。首先，训练集为以下7种类型（每个人每种类型两张图片），它们都是没有遮挡的，只是光照和表情的变化。

而测试集为以下几种类型：

由此可见，分类器的主要挑战在于遮挡。因为光照不均是出现在训练集里的，而遮挡并没有。而降低遮挡干扰的一个很重要的手段就是分块。结合以下这张图片为例可以更好地说明：

对于这两张测试集图片。戴墨镜和戴围巾显然会导致重构效果差，因为显然用字典里原子的线性组合生成一片黑色十分困难。如果分块了，虽然在围巾上和墨镜上的块同样也很难被字典表达，但是没有遮挡的部分则可以比较好得被正确表达。也就是说，假设分成9块，哪怕有5块覆盖在围巾上，4块是在未遮挡的脸上。那么因为5块围巾上的预测结果很可能都不正确，且较分散（如这5块预测的都是不同的5类，均不正确）。而未遮挡的四块其中有3块预测正确了，那么根据投票计数最高原则，最终结果将是正确的。

1. 次要原因：降采样到12*10可能过低。

5.2程序正确性验证

这是询问助教老师获得的机智的思路，即测试样本也从训练集中取，应该得到准确率为1。结果确实如此。图中，左侧列为ground truth，右侧列为9个patch投票结果计数的最大值。如(‘m-016’,8)，就意味着票数最多的类为m-016，共计8票。

因为每个样本要4s，所以并没有全部运行完，只运行了17%左右，17*12个测试样本都预测正确（虽然可能有1个patch预测失败，这可能是因为该patch覆盖着遮挡部分，一片黑色，错误是可容忍的，会被投票机制纠正）。因此可以基本断言，程序是没有错误的。

5.3尝试加速

在本次实验中，每个测试样本的分类接近4s，总共1200个测试样本就需要4800s，较为耗时。首先，我用PyCharm的profile功能，查看每个函数的执行时间占比。

其中，f是我定义的omp算法的函数，占用了95.7%的时间。pinv占用了79%的时间。（第三行的built-in method numpy.core…是包括pinv的，而pinv又是包括svd的，因为伪逆是通过SVD求得的。这两点从call graph可以看出，部分call graph如下：）。

也就是说，计算伪逆很耗时。

因此我做了加速的尝试。思路有两个，第一个是增加CPU的并行性，第二个是使用GPU。

第一个思路，没有找到可行的库，且因为numpy底层本来就有一定程度的并行，运行时的CPU占用率如下，相对较高，因此我除了修改np.dtype之外没有再尝试改进。我希望直接使用GPU，获得更大程度的加速。

（CPU实际型号为I7-8700，GPU型号为GTX 1070TI）

因此尝试使用GPU进行加速。虽然最终都没有获得更快的速度，但还是把结果记录下来：

5.3.1 性能基准：Numpy+MKL

NumPy在编译时是使用了MKL的。它是Intel提供的一个向量运算加速库。在最耗时的求伪逆这步，需要0.0001-0.006s. （因为xk是从1维增长到1400维的，所以时间增长是正常的）

而每个测试样本需要9*1400次求伪逆（同样，每次求伪逆的矩阵大小不一样）。

我试图修改了dtype为np.float32，与np.float64对比，得到如下结果：1、所有dtype=np.float64或np.float,耗时3.6s； 2、所有dtype=np.float32，耗时5.1s。

Float32居然比float64慢! 查阅资料得知，这是因为我的CPU是原生支持64位浮点数的，而32位却是软件支持（非硬件原生）。

5.3.2使用CuPy+CUDA试图加速求伪逆

CuPy是对numpy的GPU移植，试图提供与Numpy一样的API。我认为这可能会比较方便，所以先尝试了它。

代码如下：大体思路就是在每次求伪逆之前先把转移到GPU内存，然后求完伪逆之后马上转到系统内存，完成算法接下来的步骤。

每次求伪逆的时间居然增长到0.001-0.08，总共需要72s一个样本，比NumPy慢很多！我初步猜测这可能是因为将矩阵从内存到GPU的迁移会很耗时，但是哪怕我只计求伪逆的时间，不计时迁移的时间，也比CPU上操作要慢。更不要说转移是无法避免的，所以我放弃了这一条路，认为可能是CuPy库本身的问题。毕竟这是一个小众的库。

5.3.3 使用PyTorch+CUDA试图加速

上述结果令我非常震惊，我决定再使用PyTorch来试图加速。虽然它主要是被用作深度学习，但是有帖子分析到它的向量运算比Numpy要快。此外，为了避免频繁的CPU和GPU之间的通信，我把OMP稍微进行了改写，把整个OMP都放到GPU上运行，也就是说，OMP涉及到的所有矩阵和向量，包括字典，选取的列子集，求解的，残差，都尽可能早（为了减少迁移次数）地放到了GPU上，代码如下：

没想到的是，GPU占用率非常低，只有10%左右，反而CPU占用率却和NumPy+MKL差不多，整体时间也变长，从4s到20s！

这令我很疑惑，我会记住这两个GPU比CPU慢的情形，等到以后知识储备更加扎实了之后希望可以解答。我的猜测是：奇异值分解在GPU上运行的效率很低，GPU可能更加适合矩阵乘法之类的运算。

5.4 分块后的测试

因为算法本身已经固化了，所以能调整的只有降采样的尺寸以及分块的方式。

5.4.1 测试1

按照降采样到45*42，如下分块进行测试：

原图：

降采样并且分块后：

每类（人）的准确率：该类中分类正确的测试样本数/12：

总的准确率：75%（所有分类正确的测试样本数/1200）.

准确率并不是很满意，这可能是因为图片拉伸有点明显。并且特征最明显的眼睛有点被割裂。

5.4.2测试2

改变划分方式，包含完整的眼睛：

正确率按组画出如图：

总的正确率提升到了：84.4%

5.4.3测试3

之前为了方便，我都是均分成小的patch。这是为了字典可以是的形式。但实际上，如果每个patch的大小不一样，即每个sample的维度不相等，也是可以操作的。由两种办法：

1. 建立patch数个不同的字典，每个字典的为。
2. 仍然建立一个三维的字典，sample维数低的统一在末尾添0，补齐到与最大的patch的像素数相等。稍加分析可知，这是不会影响OMP算法的求解的，因为加0补齐相当于在图像上增加黑色像素。

虽然法2代码简单，但是会更为耗时，因此采用法1.构建字典如下：

我连续观察了几个人的图像，发现五官的位置在图中还是比较稳定的，这意味着可以对着五官来划分。因此，我将整个图像分为八份，例子如下：

为了节约时间，图像被进一步降采样到41*30.预测一个样本大约需要2.9s提升1s左右。

分割锚点如下：

其余代码也有少量修改，在此就不详细放出。

每类准确率图：

总体准确率：84.6%。可见，在加快速度的同时，分类准确率没有降低。

但是为什么准确率不稳定，有些人的准确率会特别低？我认为这可能和分割有关系。也就是说，准确率较低的样本，可能五官稍有位移，导致求解线性组合的时候并不准确。观察发现确实如此。

以上就是做的全部测试。因为一次测试需要一个多小时，所以就只做了这三个测试。毕竟，算法和程序框架的实现更为本质。

1. 小结

在本次课程大作业的完成过程中，我先阅读了SRC的经典论文，以及关于OMP算法的很多课件和博客。这途中，因为想对算法本质有更深入的了解，我又花了一天的时间来弥补线性代数的知识漏洞，了解基于奇异值分解的伪逆可以解决奇异矩阵无法求逆的问题，等等。也把一些不清楚的问题弄懂，加深了对算法的理解。

在实现过程中，我力求代码简单易懂不出错，充分利用了python的特性，并且使用工具（timer和Pycharm professional的profile的功能）进行了运行速度的分析和优化。也尝试将OMP算法移植至GPU，虽然最终没有得到提高，但也留下了问题，以待今后思考。

我还通过训练集自测的方法验证了代码的正确性，且尝试了不同patch划分方法对分类准确率的影响。最终，分类准确率达到了84.6%，速度也较快。如果降采样到更大分辨率，并且尝试一些更好的分割方式，准确率应该会继续提升。

注：提交的代码和该报告中的代码略有出入，是因为该报告中的代码时进行第一次实验时的代码，之后因为更换分块方式，代码有所更改。两处的代码都可以正常运行。

参考文献

[1] J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry, and Y. Ma, “Robust face recognition via sparse representation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 31, no. 2, pp. 210–227, 2009.

[2] Gilbert Strang, Introduction to Linear Algebra, Fourth Edition. 2012.