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:
- 通过监督方式实现,e.g. 监督学习里常用的Feature Selection 方法
- 通过非监督方式实现,e.g. auto encoder,HOG…
Eigenface 使用的是第二种approach
2. 关于Eigenface人面识别
Eigenface做人面识别的核心思想是:通过线性变换将raw data 做降为(encoding),然后再低纬度空间里做matching
Eigenface指: 对降为后的数据做decode 后的图片
这个过程与auto encoder无异
如果设:
I
I
I - 为输入的图片
ϕ
(
)
\phi()
ϕ() - encode过程
ϕ
′
(
)
\phi'()
ϕ′() - decode 过程
A
A
A - 压缩的数据
(
s
i
z
e
(
A
)
<
<
s
i
z
e
(
I
)
)
(size(A) << size(I))
(size(A)<<size(I))
则:
encoding:
ϕ
(
I
)
=
A
\phi(I)=A
ϕ(I)=A
decoding:
ϕ
′
(
A
)
=
e
i
g
e
n
f
a
c
e
\phi'(A) = eigenface
ϕ′(A)=eigenface
因为eigenface 是一种线性变换,encoding&decoding 可写为:
encoding:
W
I
=
A
WI = A
WI=A
decoding:
U
A
=
e
i
g
e
n
f
a
c
e
UA = eigenface
UA=eigenface
问:如何求
W
W
W和
U
U
U?
PCA(Principal Component Analysis)
3.关于PCA(in detail)
3.1. PCA直觉:
对高维数据做低维映射,使映射后的数据能够保持尽可能多的信息。这类似我们通过一个2D影子来描述一个3D人的样子。这个过程我们成为基变换
*什么是基:
定理:设 V V V为向量空间如果 r r r个向量 a 1 , a 2 , . . . a r ∈ V a_1,a_2,...a_r\in V a1,a2,...ar∈V,且满足
i) a 1 , a 2 , . . . a r a_1,a_2,...a_r a1,a2,...ar线性无关
ii) v v v中任意一个向量都可由 a 1 , a 2 , . . . a r a_1,a_2,...a_r a1,a2,...ar线性表示
那么 a 1 , a 2 , . . . a r a_1,a_2,...a_r a1,a2,...ar称之为 V V V的一个基
e.g. 比如一个二维空间(1,2),的基为(自然基)
a
1
=
[
1
0
]
a
2
=
[
0
1
]
a_1 =\left[\begin{matrix}1 \\0\end{matrix}\right] a_2 =\left[\begin{matrix}0 \\1\end{matrix}\right]
a1=[10]a2=[01]
*基变换:
定理:设向量空间V的一组基 a 1 , a 2 , . . . a r a_1,a_2,...a_r a1,a2,...ar 到另一组基 b 1 , b 2 , . . . b r b_1,b_2,...b_r b1,b2,...br 的过渡矩阵为 C C C , V V V中一个向量在这两组基下的坐标分别为 X X X和 Y Y Y, 则 X = C Y X=CY X=CY ,我们也称 X = C Y X=CY X=CY为坐标变换公式,同时也有 Y = C − 1 Y Y=C^{-1}Y Y=C−1Y
e.g.
a
=
(
2
,
3
)
T
a=(2,3)^T
a=(2,3)T 在自然基下的坐标为
(
2
,
3
)
(2,3)
(2,3),求在基
a
1
=
(
1
,
0
)
T
,
a
2
=
(
1
,
1
)
T
a_1=(1,0)^T,a_2=(1,1)^T
a1=(1,0)T,a2=(1,1)T下的坐标
解:
设坐标为
(
x
0
,
x
1
)
(x_0,x_1)
(x0,x1) 则
[
2
3
]
=
x
0
[
1
0
]
+
x
2
[
1
1
]
\left[\begin{matrix}2 \\3\end{matrix}\right] = x_0\left[\begin{matrix}1 \\0\end{matrix}\right] +x_2\left[\begin{matrix}1 \\1\end{matrix}\right]
[23]=x0[10]+x2[11]
解得:
x
0
=
−
1
x_0 =-1
x0=−1
x
2
=
3
x_2=3
x2=3
*什么样的基能够保持尽可能多的信息:
基所对应的variance越大越好。
如图:方差最大的基为
U
1
U_1
U1
存在问题:如果只是为寻找方差最大的基那么大部分的基都存在与
U
1
U_1
U1附近,及都与
U
1
U_1
U1线性相关,所以我们要求找到的基们不但方差要大同时线性无关
3.2.PCA目的:
寻找一组标准正交基,使原始数据在该基上做基变换的方差尽可能最大。
设
W
=
w
0
,
w
1
,
.
.
.
w
k
W= w_0,w_1,...w_k
W=w0,w1,...wk是我们找到基
(
k
<
<
s
i
z
e
(
I
)
)
(k<<size(I))
(k<<size(I))那么
encoding:
W
I
=
A
WI = A
WI=A
decoding:
W
T
A
=
e
i
g
e
n
f
a
c
e
≈
I
W^TA = eigenface \approx I
WTA=eigenface≈I*(标准正交基性质)
那么我们的问题就解决了XD
*为什么是“标准正交基”
定理:若 e 0 , e 1 , . . e n e_0,e_1,..e_n e0,e1,..en 是 V V V的一个标准正交基,那么 V V V中任一向量 a a a应由 e 1 , . . . , e r e_1,...,e_r e1,...,er线性表示,设表示为:
a = λ 1 e 1 + λ 2 e 2 + . . . + λ r e r a = \lambda_1e_1+\lambda_2e_2+...+\lambda_re_r a=λ1e1+λ2e2+...+λrer
为求其中的系数 λ i ( i = 1 , . . . , r ) \lambda_i(i=1,...,r) λi(i=1,...,r),有
e i T a = λ i e_i^Ta =\lambda_i eiTa=λi同样 a = e i λ i a=e_i\lambda_i a=eiλi
3.3.寻找标准正交基:
3.3.1 协方差矩阵
既然要寻找方差最大的正交基,那我们首先需要计算数据的方差variance:
V
a
r
(
a
)
=
1
m
−
1
∑
(
a
i
−
μ
)
2
Var(a)=\frac{1}{m-1}\sum(a_i-\mu)^2
Var(a)=m−11∑(ai−μ)2
同时我们需要保证两两基之前线性无关及协方差为0
C
o
v
(
a
,
b
)
=
1
m
−
1
∑
a
i
b
i
Cov(a,b)=\frac{1}{m-1}\sum a_ib_i
Cov(a,b)=m−11∑aibi
这里就引入协方差矩阵:
设
X
=
(
a
1
,
a
2
,
.
.
.
,
a
m
b
1
,
b
2
,
.
.
.
,
b
m
)
X = \left(\begin{matrix}a_1,a_2,...,a_m \\b_1,b_2,...,b_m\end{matrix}\right)
X=(a1,a2,...,amb1,b2,...,bm)
那么他的协方差矩阵为:
c
o
v
m
a
x
=
1
m
−
1
X
X
T
=
(
1
m
−
1
a
i
2
,
1
m
−
1
a
i
b
i
1
m
−
1
a
i
b
i
,
1
m
−
1
b
i
2
)
covmax = \frac{1}{m-1}XX^T=\left(\begin{matrix}\frac{1}{m-1}a_i^2,\frac{1}{m-1}a_ib_i\\\frac{1}{m-1}a_ib_i,\frac{1}{m-1}b_i^2\end{matrix}\right)
covmax=m−11XXT=(m−11ai2,m−11aibim−11aibi,m−11bi2)
这里要注意我们所用的
X
X
X的
μ
=
0
\mu=0
μ=0
可以看出协方差矩阵的对角线是方差,而其他地方是协方差,且是个对称矩阵
这里我们需要协方差都为0,我们可以通过线性变换将原协方差矩阵变为一个对角矩阵,及对角化
定理: 设A为n阶对称矩阵,则标有正交矩阵P,使 P − 1 A P = P T A P = Λ P^{-1}AP=P^TAP=\Lambda P−1AP=PTAP=Λ,其中 Λ \Lambda Λ是以A的n个特征的对角元的对角矩阵
(这里的 P P P也是个过渡矩阵,证明?)
定理:设 A A A是个 n n n阶矩阵,如果数 λ \lambda λ和 n n n维非 0 0 0列向量 x x x使关系式: A x = λ x Ax=\lambda x Ax=λx
那么, λ \lambda λ称为矩阵 A A A的特征值(eigen value), x x x称之为 A A A的对应于特征值 λ \lambda λ的特征向量(eigen vector)
由定理可知, Λ \Lambda Λ 是covmatrix的eigenvalue组成的对角矩阵。而 P P P是由eigenVector组成的矩阵。我们可以对P里的每一个vector做标准正则化,的到标准正交基
3.3.1 求标准正交基
第一步:求出对称矩阵A的eigen value 和 eigen vector
由定理可知一个矩阵
A
A
A可以写成
A
x
=
λ
x
Ax=\lambda x
Ax=λx及
(
A
−
λ
E
)
x
=
0
(A-\lambda E)x=0
(A−λE)x=0 ,
E
E
E是单位矩阵
那么它有非0解得充分必要条件是:
d
e
t
(
A
−
λ
E
)
=
0
det(A-\lambda E)=0
det(A−λE)=0
由此我们可以的到
λ
\lambda
λ继而解得x
第二部: 对eigen vector 做单位化的到
P
P
P, 使
P
T
A
P
=
P
−
1
A
P
=
Λ
P^TAP=P^{-1}AP=\Lambda
PTAP=P−1AP=Λ
设
X
=
(
x
0
,
x
1
,
.
.
x
m
)
X= (x_0,x_1,..x_m)
X=(x0,x1,..xm)为A的eigen vecors
则
P
=
(
p
0
,
p
1
,
.
.
.
,
p
m
)
=
(
1
∣
∣
x
1
∣
∣
x
1
,
1
∣
∣
x
2
∣
∣
x
2
,
.
.
1
∣
∣
x
m
∣
∣
x
m
)
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)
P=(p0,p1,...,pm)=(∣∣x1∣∣1x1,∣∣x2∣∣1x2,..∣∣xm∣∣1xm)
注:必须保证每个eigen vector都是独一无二的。
如果上述条件不成立我们需要通过施密特正交化实现
定理:设 λ i , λ j \lambda_i,\lambda_j λi,λj是对称矩阵A的两个特征值eigen values, p i , p j p_i,p_j pi,pj 是对应的特征向量eigen vectors
若 λ i ≠ λ j \lambda_i \neq \lambda_j λi̸=λj则 p i , p j p_i,p_j pi,pj正交
4. eigenface人面识别流程
reference:Eigenface for Face Detection/Recognition
算法主要分三步:
step1:训练集求出K个标准正交基
W
W
W
step2:encode训练集合数据
W
I
=
A
WI = A
WI=A
step3:encode测试集合数据 通过KNN 对测试集分类
关于eigenface
e
i
g
e
n
f
a
c
e
=
W
T
A
eigenface = W^TA
eigenface=WTA
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
I_1,I_2,..I_m
I1,I2,..Im(训练集),将原始面部图像reshape数据成
N
2
∗
1
N^2*1
N2∗1矩阵 (注 图像不一定都是方阵)
step 2 计算
m
e
a
n
f
a
c
e
=
1
m
∑
I
i
meanface = \frac{1}{m}\sum I_i
meanface=m1∑Ii 为了后面方便计算协方差矩阵
step3 数据减去 mean face:
Φ
i
=
I
i
−
m
e
a
n
f
a
c
e
\Phi_i = I_i - meanface
Φi=Ii−meanface
step4 计算协方差矩阵
C
C
C:
C
=
X
X
T
(
N
2
∗
N
2
m
a
t
r
i
x
)
X
=
[
Φ
1
,
Φ
2
,
.
.
.
,
Φ
m
]
(
N
2
∗
M
m
a
t
r
i
x
)
C=XX^T \quad (N^2*N^2 matrix) \\ X = [\Phi_1,\Phi_2,...,\Phi_m] \quad (N^2*M matrix)
C=XXT(N2∗N2matrix)X=[Φ1,Φ2,...,Φm](N2∗Mmatrix)
step5 计算标准正交基
W
W
W:
方法一 通过
W
−
1
C
W
=
W
T
C
W
=
Λ
W^{-1}CW=W^TCW=\Lambda
W−1CW=WTCW=Λ 直接计算得出正交矩阵
W
W
W和对角矩阵
Λ
\Lambda
Λ
方法二 通过
C
x
i
=
λ
i
x
i
Cx_i =\lambda_i x_i
Cxi=λixi逐个求出每个eigen value 和 vector。然后再对
x
i
x_i
xi做单位化
x
i
:
=
1
∣
∣
x
i
∣
∣
x
i
x_i := \frac{1}{||x_i||}x_i
xi:=∣∣xi∣∣1xi
step 6基于eignvalue的大小对 x i x_i xi做排序。 λ i \lambda_i λi越大说明 x i x_i xi所对应的variance越大。(因为 Λ \Lambda Λ是 C C C的相似矩阵)
step 7 选出最好的K个eigenvector 作为映射的基:
这里有两种选择方法:
- 粗暴的提取前K个最好的数据
- 基于表示信息的百分比来选择(如下伪代码)
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 S=XX^T S=XXT的大小为 10000x10000,差不多要0.8GB的内存。
这里有个小trick:我们知道对于一个M*N矩阵且M>N,改矩阵至多只有N-1个非0 eigenvalues,所以我们可以做如下计算:
设 我们的数据X(mean 为 0)为一个 MxN矩阵,其中M为数据的特征,N为数据的数量,且M>>N。
则它的协方差矩阵为
C
=
1
m
−
1
X
X
T
C=\frac{1}{m-1}XX^T
C=m−11XXT (大小为MxM)
我们目的是求出
C
C
C的eigen value 和vector :
λ
i
,
v
i
\lambda_i,v_i
λi,vi,使得
C
v
i
=
λ
i
v
i
Cv_i =\lambda_i v_i
Cvi=λivi
这里我们可以通过近似的方法求得一个nxn矩阵
C
′
C'
C′来得到
λ
i
,
v
i
\lambda_i,v_i
λi,vi,方法如下
C
′
=
X
T
X
C' =X^TX
C′=XTX 求出其 eigenvalue 和vector
γ
i
,
u
i
\gamma_i,u_i
γi,ui 得到:
C
u
i
=
γ
i
u
i
Cu_i =\gamma_iu_i
Cui=γiui
两边同时乘X得到
X
C
u
i
=
γ
i
X
u
i
XCu_i =\gamma_iXu_i
XCui=γiXui(
γ
i
\gamma_i
γi是矢量)
i.e.
X
X
T
X
u
i
=
γ
i
X
u
i
XX^TXu_i =\gamma_iXu_i
XXTXui=γiXui
C
X
u
i
=
γ
i
X
u
i
=
C
v
i
=
λ
i
v
i
CXu_i = \gamma_iXu_i=Cv_i=\lambda_iv_i
CXui=γiXui=Cvi=λivi
所以
λ
i
=
γ
i
X
u
i
=
v
i
\lambda_i =\gamma_i \quad Xu_i=v_i
λi=γiXui=vi
注此方法无法解决大数据问题。及M 和N 都很大,且N>>M
扩展如何解决大数据问题:batch 计算 researching 中…
4.2 encode训练集合数据 WI = A
这一步比较简单:对每个类别的每个图像我们都求出一个 A l A^l Al的集合就好( l l l指所在的label)
4.3 encode测试集合数据 通过KNN 对测试集分类
step 1 将测试集的每个图片NxN,reshape成Nx1,然后减去训练集的mean_face
step 2 encoder 训练集数据为
A
t
e
s
t
A_{test}
Atest
step 4 计算距离
d
=
∣
∣
A
t
e
s
t
−
A
l
∣
∣
d = ||A_{test}-A^l||
d=∣∣Atest−Al∣∣
关于距离公式的选择也是很有意思的课题。笔者在后面讨论
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