PLDA源代码分析(2)-PLDA_Verification

来源：互联网发布：aso优化推广编辑：程序博客网时间：2024/06/15 21:53

说明：此处的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相关应用

基本问题

Face verification是给出两张脸，我们判断它是否属于同一个人(不管这个人是谁)，如下图所示两个样本x1和xp，那么一共有两个模型，M0表示他们不是同一个人，M1表明他们来自同一个人(由相同identity variable h产生)。

该程序就是计算两个样本属于同一个人的概率和属于不同人概率之比，若这个概率越大，我们认为他们属于同一个人的概率也越大。

该文件有三个子程序:

PLDA_Verification. 主程序，计算这个概率之比。
preProcessPLDAModel和preProcessPLDAData是预处理程序，为了子程序getLogLikeMatchPLDA方便。
getLogLikeMatchPLDA. 计算参数X(可能只有一个样本)中的样本属于同一个人的概率。

PLDA_Verification

主程序。前半部分是调用preProcessPLDAModel和preProcessPLDAData对数据预处理，减少后面计算概率的复杂度。其中HIGHEST_N表示最多有多少个样本来自同一人。

HIGHEST_N = 2; factorModel = preProcessPLDAModel(F, G, Sigma, HIGHEST_N);Data1P = preProcessPLDAData(factorModel, Data1);Data2P = preProcessPLDAData(factorModel, Data2);

根据PLDA模型，我们知道样本 $X_{1...n}$ 属于同一个人的概率为：

$Pr(x_{1...n})=G_{x^{'}}[\mu^{'},AA^{T}+\Sigma^{'}]$

那么如上图所示，M0模型时，样本属于不同的人，那么样本在M0模型下的概率为：

$Pr(x_{1},x_{p}\mid M_{0})=Pr(x_{1}\mid M_{0})Pr(x_{p}\mid M_{0})$

那么在M1模型的情况下，样本属于同一个人，它的概率为：

$Pr(x_{1},x_{p}\mid M_{0})=Pr([x_{1},x_{2}]\mid M_{0})$

那么下面代码就是计算样本符合M1模型的概率和符合M0模型的概率之比(取对数，那么乘变为加，除变为减)

logLikeNoMatch = getLogLikeMatchPLDA(factorModel, Data1P)...    + getLogLikeMatchPLDA(factorModel, Data2P);logLikeMatch = getLogLikeMatchPLDA(factorModel, [Data1P, Data2P]);LogLikeRatio = logLikeMatch - logLikeNoMatch;

getLogLikeMatchPLDA

先看这个计算概率的程序。我们知道 $Pr(x_{1...n})=G_{x^{'}}[\mu^{'},AA^{T}+\Sigma^{'}]$ ，对概率值取对数我们可以得到：

$ln(Pr(x_{1...n}))=\frac{1}{(2\pi )^{nd}\left | AA^{T}+\Sigma^{'} \right |}e(-\frac{1}{2}x^{'T}(AA^{T}+\Sigma^{'})x^{'})$

$http://latex.codecogs.com/gif.latex?\begin{align*}%20%26%3D%20-\frac{nd}{2}(ln(2\pi))\\%20%26-\frac{1}{2}(ln(det(AA^{T}+\Sigma^{%27}))))\\%20%26-\frac{1}{2}x^{%27T}(AA^{T}+\Sigma^{%27})^{-1}x^{%27}\\%20\end{align*}$

这里x'是一个维数是nd的列向量。前面两项 $-\frac{nd}{2}(ln(2\pi))-\frac{1}{2}(ln(det(AA^{T}+\Sigma^{'})))$ 是一个常量，在程序中是factorModel.constTerm{N_DATA}.我们在子程序preProcessModel中计算。最后一项 $-\frac{1}{2}x^{'T}(AA^{T}+\Sigma^{'})^{-1}x^{'}$ 需要计算 $(AA^{T}+\Sigma^{'})^{-1}$ ，这个矩阵大小会随着样本数量而成倍增加。直接计算有难度我们对它进行如下的变换：

$(AA^{T}+\Sigma^{'})=\begin{bmatrix} FF^{'}+GG^{'} +\Sigma & FF^{'} & \cdots & FF^{'}\\ FF^{'} & FF^{'}+GG^{'}+\Sigma & \cdots & FF^{'}\\ \vdots & \vdots & \ddots & FF^{'}\\ FF^{'} & FF^{'} & \cdots & FF^{'}+GG^{'}+\Sigma \end{bmatrix}=\begin{bmatrix} FF^{'} & FF^{'} & \cdots & FF^{'}\\ FF^{'} & FF^{'} & \cdots & FF^{'}\\ \vdots & \vdots & \ddots & \vdots\\ FF^{'} & FF^{'} & \cdots & FF^{'} \end{bmatrix} + \begin{bmatrix} GG^{'} +\Sigma & 0 & \cdots & 0\\ 0 & GG^{'}+\Sigma & \cdots & 0\\ \vdots & \vdots & \ddots & 0\\ 0 & 0 & \cdots & GG^{'}+\Sigma \end{bmatrix}$

我们设变量 $\widetilde{F}$ 和 $\widetilde{\Sigma}$ 如下所示：

$\widetilde{F}=\begin{bmatrix} F\\ F\\ \vdots\\ F \end{bmatrix},\; \widetilde{\Sigma}=\begin{bmatrix} GG^{'}+\Sigma & 0 & \cdots & 0\\ 0 & GG^{'}+\Sigma & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & GG^{'}+\Sigma \end{bmatrix}$

那么根据矩阵逆定理(Binomial Inverse Theorem):

$\left(\mathbf{A}+\mathbf{UBV}\right)^{-1}=\mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{U}\left(\mathbf{B}^{-1}+\mathbf{VA}^{-1}\mathbf{U}\right)^{-1}\mathbf{VA}^{-1}.$

最终 $(AA^{T}+\Sigma^{'})^{-1}$ 可以如下计算：

$(AA^{T}+\Sigma^{'})^{-1}=(\widetilde{F}I{\widetilde{F}}'+\widetilde{\Sigma})^{-1}=\widetilde{\Sigma}^{-1}-\widetilde{\Sigma}^{-1}\widetilde{F}(I+{\widetilde{F}}'\widetilde{\Sigma}^{-1}\widetilde{F})^{-1}{\widetilde{F}}'\widetilde{\Sigma}^{-1}$

那么:

$\begin{align*} {x}'^{T}(AA^{T}+{\Sigma}')^{-1}{x}' &= \sum_{i=1}^{n}{x_{i}}'(G{G}'+\Sigma)^{-1}x_{i} \\ &+ {(\sum_{i=1}^{n}{F}'(G{G}'+\Sigma)^{-1}x_{i})}'(n{F}'(G{G}'+\Sigma)F+I)^{-1}(\sum_{i=1}^{n}{F}'(G{G}'+\Sigma)^{-1}x_{i})\end{align*}$

其中：

${x_{i}}'(G{G}'+\Sigma)^{-1}x_{i}$ 对应代码中的dataPP.quadTerm，在preProcessPLDAData子程序中计算。
$\sum_{i=1}^{n}{F}'(G{G}'+\Sigma)^{-1}x_{i}$ 对应代码中的sumWeightedData,它由dataPP(data).FTinvSx(在preProcessPLDAData子程序中计算)累加得到。
$(n{F}'(G{G}'+\Sigma)F+I)^{-1}$ 对应代码中的factorModel.invNFSPlusDiag{N_DATA}，在preProcessPLDAModel子程序中得到。

最终计算概率的代码如下：

logLike = factorModel.constTerm{N_DATA};sumWeightedData = zeros(N_HIDDEN, 1);logTerm = 0;for cData = 1 : N_DATA    logTerm = logTerm + dataPP(cData).quadTerm;            sumWeightedData = sumWeightedData + dataPP(cData).FTinvSx;endlogLike = logLike - 0.5 *(logTerm - ...    (sumWeightedData' * factorModel.invNFSFPlusIDiag{N_DATA} * sumWeightedData));return

preProcessPLDAModel

该主程序主要是对PLDA模型进行预处理计算一些变量的值如下表格：

变量名变量值invCovDiag $\Sigma^{-1}$ GWeighted $\Sigma^{-1}G$ factorModel.invTerm $({G}'\Sigma^{-1}G+I)^{-1}$ factorModel.FTranspJ ${F}'(G{G}'+\Sigma)^{-1}$ factorModel.AInv $\Sigma^{-1}$ factorModel.U $-(\Sigma)^{-1}G$ factorModel.V $({G}'\Sigma^{-1}G+I)^{-1}{G}'\Sigma^{-1}$ factorModel.MDLInvTerm $I-({G}'\Sigma^{-1}G+I)^{-1}{G}'\Sigma^{-1}G$ factorModel.logDetJ $log\left |(G{G}'+\Sigma)^{-1} \right |$ factorModel.FFfactorModel.GGfactorModel.FWeighted $(G{G}'+\Sigma)^{-1}F$ factorModel.GWeighted $\Sigma^{-1}G$ factorModel.InvNFSFPlusDiag{n} $(n{F}'(G{G}'+\Sigma)^{-1}F+I)^{-1}$ logdetInvNFSDPlusDiag $log|(n{F}'(G{G}'+\Sigma)^{-1}F+I)^{-1}|$ factorModel.constTerm{n} $-\frac{nd}{2}ln(2\pi )-\frac{1}{2}ln|AA^{T}+{\Sigma}'|$ 看几个加粗的变量：

factorModel.FTranspJ同样用的是矩阵逆定理(Binomial Inverse Theorem):

${F}'\Sigma^{-1}-{F}'\Sigma^{-1}G({G}'\Sigma^{-1}G+I)^{-1}{G}'\Sigma^{-1}={F}'(\Sigma^{-1}-\Sigma^{-1}G({G}'\Sigma^{-1}G+I)^{-1}{G}'\Sigma^{-1})={F}'(G{G}'+\Sigma)^{-1}$

代码如下：

factorModel.FTranspJ = (F.*repmat(invCovDiag,1,N_HIDDEN_DIM))'-(F'*GWeighted)*factorModel.invTerm*GWeighted';

factorModel.logDetJ用到矩阵逆定理(Binomial Inverse Theorem)和矩阵行列式定理(Matrix determinant lemma).其中矩阵行列式定理如下：

$\operatorname{det}(\mathbf{A}+\mathbf{UV}^\mathrm{T}) = \operatorname{det}(\mathbf{I} + \mathbf{V}^\mathrm{T}\mathbf{A}^{-1}\mathbf{U})\operatorname{det}(\mathbf{A}).$

那么该变量可如下计算：

$\begin{align*} log(det((G{G}'+\Sigma)^{-1})) &= log(det(I+{G}'\Sigma^{-1}G)^{-1})+log(det(\Sigma^{-1}))\\ &= log(det(I-{G}'(G{G}'+\Sigma)^{-1}{G}))+log(det(\Sigma^{-1}))\\ &= log(det(I-(I+{G}'\Sigma^{-1}G)^{-1}{G}'\Sigma^{-1}G)))+log(det(\Sigma^{-1})) \end{align*}$

其中最后一步利用公式： $(P^{-1}+{H}'R^{-1}H)^{-1}{H}'R^{-1}=P{H}'(HP{H}'+R)^{-1}$ 。

具体代码如下所示：

[U L V]=svd(factorModel.MDLInvTerm);factorModel.logDetJ =sum(log(diag(L))) + sum(log(invCovDiag));

factorModel.InvNFSFPlusDiag{n}参见SVD和PCA
factorModel.constTerm{n}后半部分 $ln(det(A{A}'+{\Sigma}'))$ ,可以如下计算：

$\begin{align*} ln(det(A{A}'+{\Sigma}')) &= ln(det(\begin{bmatrix} F{F}'+G{G}'+\Sigma & F{F}' & \cdots & F{F}'\\ F{F}' & F{F}'+G{G}'+\Sigma & \cdots & F{F}'\\ \vdots & \vdots & \ddots & \vdots\\ F{F}' & F{F}' & \cdots & F{F}'+G{G}'+\Sigma \end{bmatrix}))\\ &= ln(det(\begin{bmatrix} nF{F}'+G{G}'+\Sigma & F{F}' & \cdots & F{F}'\\ 0 & G{G}'+\Sigma & \cdots & F{F}'\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & G{G}'+\Sigma \end{bmatrix}))\\ &= ln(det(nF{F}'+G{G}'+\Sigma)) + (n-1)ln(det(G{G}'+\Sigma))\\ &= -ln(det(n{F}'(G{G}'+\Sigma)^{-1}F+I)^{-1})-ln(det(G{G}'+\Sigma)^{-1}) - (n-1)ln(det(G{G}'+\Sigma)^{-1})\\ &= -ln(det(n{F}'(G{G}'+\Sigma)^{-1}F+I)^{-1})-nln(det(G{G}'+\Sigma)^{-1}) \end{align*}$

分别对应于logdetInvNFSDPlusDiag和factorModel.logDetJ，代码部分如下所示：

for cN = 1 : HIGHEST_N    [U, L, V] = svd(cN*F'*factorModel.FWeighted+eye(N_HIDDEN_DIM));    DiagL = diag(L);    factorModel.invNFSFPlusIDiag{cN} = V * diag(1 ./ DiagL) * U';        logdetInvNFSFPlusIDiag = sum(log(1 ./ DiagL));            factorModel.constTerm{cN} = - (cN * N_DATA_DIM / 2) * log(2* pi)...        + cN / 2 * factorModel.logDetJ + 0.5 * logdetInvNFSFPlusIDiag;end

preProcessPLDAData

该子程序主要计算 ${x_{i}}'(G{G}'+\Sigma)^{-1}x_{i}$ ，对应程序中的quadTerm。根据矩阵逆定理(Binomial Inverse Theorem)可以如下计算：

${x_{i}}'(G{G}'+\Sigma)^{-1}x_{i}={x_{i}}'\Sigma^{-1}x_{i}-{x_{i}}'\Sigma^{-1}G({G}'\Sigma^{-1}G+I)^{-1}{G}'\Sigma^{-1}x_{i}$

分别对应程序中的quadTerm1和quadTerm2，代码如下：

for cData = 1 : N_DATA    quadTerm1 = (factorModel.AInv.*data(:,cData))'*data(:,cData);    quadTerm2 = (data(:,cData)'*factorModel.GWeighted)*factorModel.invTerm*(factorModel.GWeighted'*data(:,cData));    quadTerm = quadTerm1-quadTerm2;    dataPP(cData).quadTerm = quadTerm;    end