PLDA源代码分析(1)-PLDA_Train

来源：互联网发布：log4cpp linux编译编辑：程序博客网时间：2024/05/18 01:23

本博文转载于http://write.blog.csdn.net/postedit?ref=toolbar

说明：此处的LDA对应于Linear Discriminant Analysis，PLDA即对应于Probabilistic LDA. 该代码对应的文章为ICCV2007 paper Probabilistic Linear Discriminant Analysis for Inferences About Identity，源代码可以从 Prince Vision Lab处下载。虽然源码虽然不长结构比较清楚，但是运用到了一定的矩阵知识，所以对源码分析稍作分析。

1、PLDA 训练(Training)源码分析

2、PLDA识别(Recognition)源码分析

3、PLDA相关应用

基本问题

PLDA的基本模型如下：

$Pr(x_{i,j}|h_{i},w_{i,j},\theta )=G_{x}(\mu +Fh_{i}+Gw_{i,j},\Sigma )$

$Pr(h_{i})=G_{h}(0,I)$

$Pr(w_{i,j})=G_{h}(0,I)$

Learning LIV Models给出数据集 $x_{i,j}$ ，利用EM算法求得参数 $\theta=\{\mu,F,G,\Sigma \}$ ：

E-Step:计算两个期望：

$E[y_{i}]=(A^{T}\Sigma^{'-1}A+I)^{-1}A^{T}\Sigma^{'-1}(x_{i}-\mu^{'})$

$E[y_{i}y_{i}^{T}]=(A^{T}\Sigma^{'-1}A+I)^{-1}+E[y_{i}]E[y_{i}]^{T}$

其中：

$A=\begin{bmatrix} F & G & 0 & \cdots & 0\\ F & 0 & G & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ F & 0 & 0 & \cdots & G\\ \end{bmatrix}$

$\Sigma=\begin{bmatrix} \Sigma & \cdots & 0\\ \vdots & \ddots & \vdots\\ 0 & \cdots & \Sigma\\ \end{bmatrix}$

M-Step:计算更新三个参数：

$\mu=\frac{1}{IJ}\sum_{i,j}x_{i,j}$

$B=(\sum_{i,j}(x_{i,j}-\mu)E[z_{i}]^{T})(\sum_{i,j}E[z_{i}z_{i}]^{T})^{-1}$

$\Sigma=\frac{1}{IJ}\sum_{i,j}Diag[(x_{i,j}-\mu)x_{i,j}-\mu)^{T}-BE[z_{i}](x_{i,j}-\mu)^{T})]$

其中均值不变，Diag操作是取矩阵对角元素构成一个对角阵，B和 $z_{i}$ 分别如下所示：

$z_{i,j}=\begin{bmatrix} h_{i}\\ w_{i,j}\\ \end{bmatrix}$ ， $B=[F,G]$

PLDA_Train.m共有三个函数，其中PLDA_Train是EM算法的主函数，getExpectedValuePLDA是计算期望的子函数，trainPCA是初始化F的PCA过程。

PLDA_Train

输入参数：

Data:d*n的训练数据，d为样本维数，n为样本个数。
ImageID:n*m的一个稀疏矩阵，m表示有多少个人，若ImageID[i][j]=1那么表示第i个样本来自第j个人。
N_ITER:EM算法迭代次数
N_F，N_G:因子个数，即矩阵F和G的第二维(其第一维是d)

EM算法迭代过程：

[plain] view plaincopy
for cIter = 1 : N_ITER  
    %E-Step  
    %M-Step  
end  

更新 $E[y_{i}]$ 和 $E[y_{i}y_{i}^{T}]$ ，这里Eh是 $[E[y_{1}],E[y_{1}],...,E[y_{n}]]$ 和EhhSum是 $\sum_{i,j}E[z_{i}z_{i}^{T}]$ ，二者均为矩阵。

[plain] view plaincopy
[Eh EhhSum] = getExpectedValuesPLDA(F, G, Sigma, Data, ImageID);  

更新B

[plain] view plaincopy
xhSum = zeros(size(Data,1), size(Eh, 1));  
for cData = 1 : N_DATA  
   xhSum = xhSum + Data(:, cData) * Eh(:, cData)';  
end;  
   FGEst = xhSum * inv(EhhSum) ;  

更新 $\Sigma$

若Data如下所示：

$\begin{bmatrix} x_{1,1} & x_{2,1} & \cdots & x_{n,1}\\ x_{1,2} & x_{2,2} & \cdots & x_{n,2}\\ \vdots & \vdots & \ddots & \vdots\\ x_{1,d} & x_{2,d} & \cdots & x_{n,d} \end{bmatrix}$

那么mean(Data*Data',2)如下所示：

$\begin{bmatrix} \frac{1}{n}(x_{1,1}^{2} + x_{2,1}^{2} + \cdots + x_{n,1}^{2})\\ \frac{1}{n}(x_{1,2}^{2} + x_{2,2}^{2} + \cdots + x_{n,2}^{2})\\ \vdots \\ \frac{1}{n}(x_{1,d}^{2} + x_{2,d}^{2} + \cdots + x_{n,d}^{2}) \end{bmatrix}=\frac{1}{IJ}\sum_{i,j}Diag(x_{i,j}-\mu)(x_{i,j}-\mu)^{T}$

后面一部分也类似。

[plain] view plaincopy
Sigma = mean(Data .* Data - (FGEst * Eh) .* Data,2);  

根据B的定义取得F和G

[plain] view plaincopy
F = FGEst(:, 1 : N_F);  
G = FGEst(:, N_F + 1 : end);  

getExpectedValuesPLDA

主要思想：

该子程序主要计算 $[E[y_{1}],E[y_{1}],...,E[y_{n}]]$ 和 $\sum_{i,j}E[z_{i}z_{i}^{T}]$ ，二者均为矩阵。前半部分预先计算 $(A^{T}\Sigma^{'-1}A+I)^{-1}$ 即样本个数为repeatValues(repeatValues为每个人的样本数目)的指，后半部分具体计算这两个矩阵。

预处理打表

计算表的大小，如PLDA_Train_Demo.m中repeatValues=[2,3,4]共三种情况，那么nRepeatValues=3为循环的个数，invTermsAll存储 $(A^{T}\Sigma^{'-1}A+I)^{-1}$ 的值。

[plain] view plaincopy
repeatValues = unique(sum(ImageID));  
nRepeatValues = length(repeatValues);  
invTermsAll = cell(nRepeatValues, 1);<pre name="code" class="plain"></pre>  

下面的循环计算不同repeatValues时 $(A^{T}\Sigma^{'-1}A+I)^{-1}$ 的值。

[plain] view plaincopy
for cRepeatVal = 1 : nRepeatValues  
    thisRepVal = repeatValues(cRepeatVal);  
    ...  
    invTermsAll{repeatValues(cRepeatVal)} = invTerm;  
end  

对于每个thisRepVal，易知 $A^{T}\Sigma^{'-1}A$ 是一个 (N_HID_DIM+thisRepVal*N_HID_DIM_NOISE)* (N_HID_DIM+thisRepVal*N_HID_DIM_NOISE)大小的矩阵如下：

$A^{T}\Sigma^{'-1}A=\begin{bmatrix} thisRepVal*F^{'}\Sigma F & F^{'}\Sigma G & \cdots & F^{'}\Sigma G\\ G^{'}\Sigma F & G^{'}\Sigma G & \cdots & F^{'}\Sigma G\\ \vdots & \vdots & \ddots & \vdots\\ G^{'}\Sigma F & G^{'}\Sigma G & \cdots & G^{'}\Sigma G \end{bmatrix}$

那些下面代码段就是分块对矩阵进行赋值，其中左上角矩阵维数是H_HID_DIM*H_HID_DIM，其他矩阵维数是N_HID_DIM_NOISE*N_HID_DIM_NOISE。

[plain] view plaincopy
ATISigA(1:N_HID_DIM,1:N_HID_DIM) = thisRepVal*weightedF'*F;  
for cMat = 1:thisRepVal  
    ATISigA(N_HID_DIM+(cMat-1)*N_HID_DIM_NOISE+1:N_HID_DIM+cMat*N_HID_DIM_NOISE,1:N_HID_DIM) = weightedG'*F;  
    ATISigA(1:N_HID_DIM,N_HID_DIM+(cMat-1)*N_HID_DIM_NOISE+1:N_HID_DIM+cMat*N_HID_DIM_NOISE) = weightedF'*G;  
    ATISigA(N_HID_DIM+(cMat-1)*N_HID_DIM_NOISE+1:N_HID_DIM+cMat*N_HID_DIM_NOISE,...  
        N_HID_DIM+(cMat-1)*N_HID_DIM_NOISE+1:N_HID_DIM+cMat*N_HID_DIM_NOISE) = weightedG'*G;  
end;  

有了 $A^{T}\Sigma^{'-1}A$ 剩余代码就是计算并存储 $(A^{T}\Sigma^{'-1}A+I)^{-1}$

计算期望值

外循环是每个人，内循环计算特定人所有样本的期望值，更新期望和。

[plain] view plaincopy
for cInd = 1 : N_INDIV  
    for cFaces = 1 : nFaces  
        Eh(:,thisImIndex(cFaces)) = thisEh(thisIndex);  
        EhhSum=EhhSum+thisEhh(thisIndex,thisIndex);  
    end  
end  

对每个人，首先获得样本数量以及样本所对应的位置。

[plain] view plaincopy
nFaces = full(sum(ImageID(:,cInd)));  
thisImIndex = find(ImageID(:,cInd));  

下面就是计算 $A^{T}\Sigma^{'-1}X$ ，计算的过程和计算 $A^{T}\Sigma^{'-1}A$ 类似，也是对矩阵进行分块复制。这里的X指该属于某个人的所有样本。

[plain] view plaincopy
dataAll = x(:,thisImIndex).*repmat(1./Sigma,1,nFaces);  
ATISigX = zeros(N_HID_DIM+nFaces*N_HID_DIM_NOISE,1);  
ATISigX(1:N_HID_DIM,:) = sum(F'*dataAll,2);  
for cIm = 1 : nFaces  
    ATISigX(N_HID_DIM+(cIm-1)*N_HID_DIM_NOISE+1:N_HID_DIM+cIm*N_HID_DIM_NOISE,:) = G'*dataAll(:,cIm);  
end;  

下面计算更新 $E[y_{i}]$ (thisEh)和 $E[y_{i}y_{i}^{T}]$ (thisEhh)。

[plain] view plaincopy
thisEh = invTerm*ATISigX;  
thisEhh = invTerm+thisEh*thisEh';  

最后就是for循环更新 $[E[y_{1}],E[y_{1}],...,E[y_{n}]]$ 和 $\sum_{i,j}E[z_{i}z_{i}^{T}]$ ，前者分块赋值，后者累加求和。

[plain] view plaincopy
for cFaces = 1 : nFaces  
    thisIndex = [1:N_HID_DIM  N_HID_DIM+(cFaces-1)*N_HID_DIM_NOISE+1:N_HID_DIM+cFaces*N_HID_DIM_NOISE];  
    Eh(:,thisImIndex(cFaces)) = thisEh(thisIndex);  
    EhhSum=EhhSum+thisEhh(thisIndex,thisIndex);  
end  

trainPCA

PCA过程，详见 SVD和PCA

1 0