从AutoEncoder角度理解EigenFace人面识别

1. 图片识别直觉

人面识别包括所有的图像识别，其核心是一种matching problem。及图像上的一些特征是否与我们记忆里的一些特征相匹配。比如警察在监控数据里寻找嫌疑犯的时候，往往会通过嫌疑犯的一些体貌特征来对数据进行人工匹配，图像识别也是如此。

什么样的特征来进行比较显得极为重要。一种粗暴的方法是用全部像素做pixel to pixel比较（raw data）。（计算 test 集的每一个 pixel 与 对应的training集的每一个pixel的距离。）这样如果该图片与哪一个(或者几个)个traning集图片距离最短，则认为该图片与那个training集的图片属于一个类别。这种分类算法称为 k nearest neighbors (KNN)

像素之前做matching的不足：
1， 像素太多，意味着数据的维度太大，现有的距离公式就很难有效的测量两个图片之间的距离
2， 信息太多计算机内存也会出现问题

解决方法：
对像素做特征提取/压缩/降为

笔者认为，特征提取也好，数据压缩也好，降为也好实质上指的是一种东西： 用小数据来描述大数据‘raw data’。
实现的方法：往往是让raw data 通过某些数学变换实现的。

对数据做特征提取/压缩/降为，常用的有两种approach：

1. 通过监督方式实现，e.g. 监督学习里常用的Feature Selection 方法
2. 通过非监督方式实现，e.g. auto encoder,HOG…

Eigenface 使用的是第二种approach

2. 关于Eigenface人面识别

Eigenface做人面识别的核心思想是：通过线性变换将raw data 做降为（encoding），然后再低纬度空间里做matching

Eigenface指： 对降为后的数据做decode 后的图片

这个过程与auto encoder无异

如果设：
$I$ - 为输入的图片
$\varphi ()$ - encode过程
${\varphi }^{\prime }()$ - decode 过程
$A$ - 压缩的数据 $(size(A)<<size(I))$
则：
encoding: $\varphi (I)=A$
decoding: ${\varphi }^{\prime }(A)=eigenface$
因为eigenface 是一种线性变换，encoding&decoding 可写为：
encoding: $WI=A$
decoding: $UA=eigenface$

问：如何求$W$和$U$？
PCA（Principal Component Analysis）

3.关于PCA（in detail）

3.1. PCA直觉：

对高维数据做低维映射，使映射后的数据能够保持尽可能多的信息。这类似我们通过一个2D影子来描述一个3D人的样子。这个过程我们成为基变换

*什么是基：

定理：设$V$为向量空间如果$r$个向量$a_1,a_2,...a_r\in V$,且满足
i) $a_1，a_2,...a_r$线性无关
ii)$v$中任意一个向量都可由$a_1，a_2,...a_r$线性表示
那么$a_1，a_2,...a_r$称之为$V$的一个基

e.g. 比如一个二维空间(1,2),的基为（自然基）
$$a_1=\left[\begin{array}{c}1\\ 0\end{array}\right]a_2=\left[\begin{array}{c}0\\ 1\end{array}\right]$$

*基变换：

定理：设向量空间V的一组基$a_1，a_2,...a_r$ 到另一组基$b_1，b_2,...b_r$ 的过渡矩阵为 $C$ ，$V$中一个向量在这两组基下的坐标分别为$X$和$Y$, 则$X=CY$ ,我们也称$X=CY$为坐标变换公式，同时也有$Y=C^{-1}Y$

e.g. $a=(2,3{)}^T$ 在自然基下的坐标为$(2,3)$,求在基$a_1=(1,0{)}^T,a_2=(1,1{)}^T$下的坐标
解：
设坐标为 $(x_0,x_1)$ 则
$$\left[\begin{array}{c}2\\ 3\end{array}\right]=x_0\left[\begin{array}{c}1\\ 0\end{array}\right]+x_2\left[\begin{array}{c}1\\ 1\end{array}\right]$$
解得：$x_0=-1$ $x_2=3$

*什么样的基能够保持尽可能多的信息:

基所对应的variance越大越好。
如图：方差最大的基为$U_1$

存在问题：如果只是为寻找方差最大的基那么大部分的基都存在与$U_1$附近,及都与$U_1$线性相关,所以我们要求找到的基们不但方差要大同时线性无关

3.2.PCA目的：

寻找一组标准正交基，使原始数据在该基上做基变换的方差尽可能最大。

设$W=w_0,w_1,...w_k$是我们找到基$(k<<size(I))$那么
encoding: $WI=A$
decoding: $W^{T}A=eigenface\approx I$*（标准正交基性质）
那么我们的问题就解决了XD

*为什么是“标准正交基”

定理：若$e_0,e_1,..e_n$ 是$V$的一个标准正交基，那么$V$中任一向量$a$应由$e_1,...,e_r$线性表示，设表示为：
$a={\lambda }_1{e}_1+{\lambda }_2{e}_2+...+{\lambda }_r{e}_r$
为求其中的系数${\lambda }_i(i=1,...,r)$,有
$e_i^{T}a={\lambda }_i$同样 $a=e_i{\lambda }_i$

3.3.寻找标准正交基:

3.3.1 协方差矩阵

既然要寻找方差最大的正交基，那我们首先需要计算数据的方差variance：
$$Var(a)=\frac{1}{m-1}\sum (a_i-\mu {)}^2$$

同时我们需要保证两两基之前线性无关及协方差为0
$$Cov(a,b)=\frac{1}{m-1}\sum a_i{b}_i$$

这里就引入协方差矩阵：
设 $$X=\left(\begin{array}{c}{a}_1,a_2,...,a_m\\ b_1,b_2,...,b_m\end{array}\right)$$

那么他的协方差矩阵为：
$$covmax=\frac{1}{m-1}X{X}^T=\left(\begin{array}{c}\frac{1}{m-1}{a}_i^2,\frac{1}{m-1}{a}_i{b}_i\\ \frac{1}{m-1}{a}_i{b}_i,\frac{1}{m-1}{b}_i^2\end{array}\right)$$
这里要注意我们所用的$X$的$\mu =0$

可以看出协方差矩阵的对角线是方差，而其他地方是协方差，且是个对称矩阵

这里我们需要协方差都为0，我们可以通过线性变换将原协方差矩阵变为一个对角矩阵，及对角化

定理： 设A为n阶对称矩阵，则标有正交矩阵P，使$P^{-1}AP=P^{T}AP=\Lambda$,其中$\Lambda$是以A的n个特征的对角元的对角矩阵

（这里的$P$也是个过渡矩阵，证明？）

定理：设$A$是个$n$阶矩阵，如果数$\lambda$和$n$维非$0$列向量$x$使关系式： $Ax=\lambda x$
那么，$\lambda$称为矩阵$A$的特征值（eigen value）,$x$称之为$A$的对应于特征值$\lambda$的特征向量（eigen vector）

由定理可知，$\Lambda$ 是covmatrix的eigenvalue组成的对角矩阵。而$P$是由eigenVector组成的矩阵。我们可以对P里的每一个vector做标准正则化，的到标准正交基

3.3.1 求标准正交基

第一步：求出对称矩阵A的eigen value 和 eigen vector
由定理可知一个矩阵$A$可以写成
$Ax=\lambda x$及
$(A-\lambda E)x=0$ ，$E$是单位矩阵
那么它有非0解得充分必要条件是:
$det(A-\lambda E)=0$
由此我们可以的到$\lambda$继而解得x

第二部: 对eigen vector 做单位化的到$P$, 使$P^{T}AP=P^{-1}AP=\Lambda$
设 $X=(x_0,x_1,..x_m)$为A的eigen vecors
则 $P=(p_0,p_1,...,p_m)=(\frac{1}{||x_1||}{x}_1,\frac{1}{||x_2||}{x}_2,..\frac{1}{||x_m||}{x}_m)$
注：必须保证每个eigen vector都是独一无二的。
如果上述条件不成立我们需要通过施密特正交化实现

定理：设${\lambda }_i,{\lambda }_j$是对称矩阵A的两个特征值eigen values，$p_i,p_j$ 是对应的特征向量eigen vectors
若 ${\lambda }_i̸={\lambda }_j$则$p_i,p_j$正交

4. eigenface人面识别流程

reference:Eigenface for Face Detection/Recognition

算法主要分三步：
step1：训练集求出K个标准正交基 $W$
step2：encode训练集合数据 $WI=A$
step3：encode测试集合数据 通过KNN 对测试集分类

关于eigenface
$$eigenface=W^{T}A$$

4.1 训练集求出K个标准正交基W

step0 预处理，eigenface算法要求面部图像必须是对其的，及嘴巴和眼睛的在所有数据集上位置必须保持一致。关于如何align & resize faces 可以参考Facial Landmark Detection.
也可以使用已经处理好的数据：
AT&T Facedatabase
Yale Facedatabase A
Extended Yale Facedatabase B

step1 提取数据$I_1,I_2,..I_m$(训练集)，将原始面部图像reshape数据成 $N^2\ast 1$矩阵 (注 图像不一定都是方阵)

step 2 计算$meanface=\frac{1}{m}\sum I_i$ 为了后面方便计算协方差矩阵

step3 数据减去 mean face:
$${\Phi }_i=I_i-meanface$$

step4 计算协方差矩阵$C$:
$$\begin{gathered} C=X{X}^T(N^2\ast N^{2}matrix) \\ X=[{\Phi }_1,{\Phi }_2,...,{\Phi }_m](N^2\ast Mmatrix) \end{gathered}$$

step5 计算标准正交基$W$:
方法一 通过 $W^{-1}CW=W^{T}CW=\Lambda$ 直接计算得出正交矩阵$W$和对角矩阵$\Lambda$
方法二 通过$C{x}_i={\lambda }_i{x}_i$逐个求出每个eigen value 和 vector。然后再对$x_i$做单位化
$$x_i:=\frac{1}{||x_i||}{x}_i$$

step 6基于eignvalue的大小对$x_i$做排序。${\lambda }_i$越大说明$x_i$所对应的variance越大。（因为$\Lambda$是$C$的相似矩阵）

step 7 选出最好的K个eigenvector 作为映射的基:
这里有两种选择方法:

1. 粗暴的提取前K个最好的数据
2. 基于表示信息的百分比来选择（如下伪代码）

keepInfo = 0.9 //keep 90% informations
eig_vectors = sorted_decrease(eig_vectors,by=eig_values)
total_info = sum(eig_values)
info_kept = sorted_decrease(eig_values)/total_info*100
cum_info_kept = cumsum(info_kept)

W = []
for i = 0:size(cum_info_kept) do
	wi = eig_vector[i]
	percentage = cum_info_kept[i]
	W.append[wi]
	if percentage > keepInfo*100 do
		break
	end if
end for

PCA致命问题-内存消耗巨大

从step 6 可以看出 在计算eigenvalue 和 eigenvetor的时候我们需要内存去存储：C（协方差矩阵），W（正交矩阵以及） $\Lambda$（对角矩阵）这里引用 维基百科的例子：

假设我们有400张像素为100x100的图片,那么对应的协方差矩阵$S=X{X}^T$的大小为 10000x10000,差不多要0.8GB的内存。

这里有个小trick：我们知道对于一个M*N矩阵且M>N,改矩阵至多只有N-1个非0 eigenvalues,所以我们可以做如下计算：
设 我们的数据X(mean 为 0)为一个 MxN矩阵，其中M为数据的特征，N为数据的数量，且M>>N。
则它的协方差矩阵为$C=\frac{1}{m-1}X{X}^T$ (大小为MxM)
我们目的是求出$C$的eigen value 和vector ：${\lambda }_i,v_i$,使得 $C{v}_i={\lambda }_i{v}_i$
这里我们可以通过近似的方法求得一个nxn矩阵$C^{\prime }$来得到${\lambda }_i,v_i$,方法如下
$C^{\prime }=X^{T}X$ 求出其 eigenvalue 和vector ${\gamma }_i,u_i$ 得到：
$C{u}_i={\gamma }_i{u}_i$
两边同时乘X得到
$XC{u}_i={\gamma }_{i}X{u}_i$(${\gamma }_i$是矢量)
i.e. $X{X}^{T}X{u}_i={\gamma }_{i}X{u}_i$
$CX{u}_i={\gamma }_{i}X{u}_i=C{v}_i={\lambda }_i{v}_i$
所以${\lambda }_i={\gamma }_{i}X{u}_i=v_i$

注此方法无法解决大数据问题。及M 和N 都很大，且N>>M

扩展如何解决大数据问题：batch 计算 researching 中…

4.2 encode训练集合数据 WI = A

这一步比较简单：对每个类别的每个图像我们都求出一个$A^l$的集合就好($l$指所在的label)

4.3 encode测试集合数据 通过KNN 对测试集分类

step 1 将测试集的每个图片NxN,reshape成Nx1,然后减去训练集的mean_face
step 2 encoder 训练集数据为$A_{test}$
step 4 计算距离 $d=||A_{test}-A^l||$
关于距离公式的选择也是很有意思的课题。笔者在后面讨论
step 5 分类：
方法一：我们可以手动设置threshold，及如果distence 低于某个threshold 则为类别0…
方法二：取argmax
方法三：KNN

4.4 Eigenface python 代码

#-- step 1, get data
    images,T = uti.getATface()
    train_org,trT,test_org,teT = uti.getTrainTestData(images,T) # here i use AT_face
    image_shape = images[0].shape
   	  #-- step 2, format data to NN*1 matrix
    train = toNN1Format(train_org)
    test = toNN1Format(test_org)

    #-- step 3, get mean face
    mean_face = np.mean(train,axis=1,keepdims=1)

    #-- step 4, substract the mean face
    demean_train  = train - mean_face

    #-- step 5, comput the convariance matrix C
    data_Size = train.shape
    #--- apply memory optimization process --- #
    optOn= False
    if np.argmax(data_Size) == 0: # case feature size is greater than datassize:
        print '**memory base optimization active'
        optOn = True

    if optOn: # C'= A'A
        C = np.dot(np.transpose(demean_train),demean_train)
    else: # C = AA'
        C = np.dot(demean_train,np.transpose(demean_train))

    #-- step 6, comput the eigenvector C
    if optOn: # u = Au_
        V, U_ = np.linalg.eig(C)
        U = np.dot(demean_train,U_)
    else:
        V, U = np.linalg.eig(C)
    #--- normalize ui so ||ui|| = 1
    norm_U = np.linalg.norm(U,axis=0,keepdims=1)
    U = U/norm_U

    #-- step 7,keep only K eigenvectors base on K largest eigenvalue
    #--- sort eigenvector
    eig_pairs = [(np.abs(V[i]),U[:,i]) for i in range(len(V))]
    eig_pairs.sort(key=lambda  x: x[0],reverse=True)
    #--- ref :https://towardsdatascience.com/principal-component-analysis-your-tutorial-and-code-9719d3d3f376
    tot = sum(V)
    var_exp =  [(i / tot) * 100 for i in sorted(V, reverse=True)]
    cum_var_exp = np.cumsum(var_exp)

    Uk = []
    firstRun = True
    for index, percentage in enumerate(cum_var_exp):
        ui = np.transpose(np.array([eig_pairs[index][1]]))
        if firstRun:
            Uk = ui
            firstRun = False
        else:
            Uk = np.hstack([Uk,ui])
        if percentage > keepData*100:
            break
 	#----------------------------------------------#
    #      face Representing onto this basis       #
    #----------------------------------------------#
    #-- randomly pick a face image from traning set
    #random_index = random.randrange(len(train_org))
    #a_image = train[:,random_index:random_index+1]

    random_index = random.randrange(len(test_org))
    a_image = test[:,random_index:random_index+1]

    #-- substarct mean
    demean_image = a_image - mean_face

    #-- get W , W = U'A
    W = np.dot(np.transpose(Uk),demean_image)

    #-- eigenface = WUk
    eigen_face = np.dot(Uk,W)+mean_face
    eigen_face = np.uint8(toOrgFormat(eigen_face,image_shape))
    a_image = toOrgFormat(a_image,image_shape)

    face = np.hstack([a_image,eigen_face])
    cv2.imshow('org',face)
    cv2.waitKey()
     #------------------------------------------------#
    #       face Recognition Using EigenFace         #
    #------------------------------------------------#
    #-- step 1 get set of training W crossponding with its label, Wl_list
    K = 1

    n_of_class = np.max(trT)
    Wl_list = []

    for T in xrange(n_of_class):
        l_indexes = np.nonzero(trT==T+1)[0]
        Wl = []
        image = demean_train[:,l_indexes]
        Wl = getW(image,Uk)
        Wl_list.append(Wl)
    #-- step 2 get subtract mean
    demean_test = test - mean_face

    #-- step 3 comput distence between test and train and use argmin for classification
    K = max(1,K)
    Y = []
    for i in xrange(len(test_org)):    # <- change
        W = getW(demean_test[:,i],Uk)  # <- change
        dist_list = []
        class_list =[]
        for index,Wl in enumerate(Wl_list):
            dist = getDist_list(W,Wl)#np.linalg.norm(Wl-W,axis=0,keepdims=1)
            dist = np.sort(dist)[0:K] # save memory space
            dist_list.append(dist)
            class_list.append(np.ones(len(dist))*index)

        class_list = np.hstack(class_list)
        #-- get Nearst N --#
        class_list = class_list[dist_list.argsort()[0:K]]
        if K > 1:
            label = countFrequency(class_list)+1
        else:
            label = np.argmin(dist_list)+1
        Y.append(label)

    #-- step 4 calculate accuracy
    correct = 0.0
    for index,l in enumerate(teT):   # <- change
        if l == Y[index]:
            correct += 1.0
    accuracy = correct/float(len(teT))  # <- change

    print 'correct %f'%correct
    print 'accuracy %f'%accuracy

5. AutoEncoder 人面识别

6. PCA/AutoEncoer data retrieval

7. 关于距离公式