I-Vector算法

以上的推导来源于leon晋：https://www.zhihu.com/question/63978977

变量维度声明

训练集有S句话，有$S = \{1,2，\dots,S\}$

第$s$句话有$H(s)$帧，每一帧的MFCC特征为$D$维，则每句语音的声学特征矩阵为：$X(s)=[x_1,x_2,\dots,x_i,\dots,x_{H(s)}]$

MFCC特征维度为：$(n,D)$

第k个高斯分布的参数为：$\lambda = \{\pi_k,\mu_k,\Sigma_k\}_{k=1}^K$

高斯均值超向量(列向量) $m=[\mu_1;\mu_2;\dots;\mu_c]$ ，维度大小为：($K\cdot D,1$)

综上，可以引出i-vector重要公式：

$M(s) = m + T w(s)$

1. M(s) 是UBM根据说话人MAP后面的均值超向量，之前直接用来识别了，现在也需要垒起来。
2. m 是均值超向量
3. T 是 映射空间，我把它叫做映射规则或者变换规则
4. w(s) 是$i-vector$

B-W的3种充分统计量

这部分是基于GMM高斯混合模型推导而出的B-W统计量，根据每个分布中的后验概率公式

这里计算B-W 统计量的目的是为了后面公式的简洁表达

$\begin{aligned} &N_{k}(s)=\sum_{i=1}^{H(s)} P_{\lambda}\left(k=z_i | x_{i}\right)\\ &F_{k}(s)=\sum_{i=1}^{H(s)} P_{\lambda}\left(k=z_i | x_{i}\right)\left(x_{i}-\mu_{k}\right)\\ &S_{k}(s)=\sum_{i=1}^{H(s)} P_{\lambda}\left(k=z_i | x_{i}\right)\left(x_{i}-\mu_{k}\right)\left(x_{i}-\mu_{k}\right)^{T} \end{aligned}$

注：这和EM算法迭代公式不同，主要区别在$F_k(s)$，这里做了一个去均值操作。此外，公式中需要对$F_k(s)$和$S_k(s)$除以当前$k$分布下的分布个数，即$N_k(s)$。这里为了向量化处理公式，把$N_k(s)$和$F(s)$做了矩阵扩充。

$N(s)=\left[\begin{array}{ccc} {N_{1}(s) I} & {} & {0} \\ {} & {\dots} & {} \\ {0} & {} & {N_{c}(s) I} \end{array}\right]$

$F(s)=\left[\begin{array}{c} {F_{1}(s)} \\ {\dots} \\ {F_{c}(s)} \end{array}\right]$

% 普通的充分统计量估计
[N, F] = expectation(data, ubm); % N的维度:(nmix,1)

function [N, F] = expectation(data, gmm)
% compute the sufficient statistics
post = postprob(data, gmm.mu, gmm.sigma, gmm.w(:));
N = sum(post, 2);
F = data * post';

function [post, llk] = postprob(data, mu, sigma, w)
% compute the posterior probability of mixtures for each frame
post = lgmmprob(data, mu, sigma, w);
llk  = logsumexp(post, 1);
post = exp(bsxfun(@minus, post, llk));

function logprob = lgmmprob(data, mu, sigma, w)
% compute the log probability of observations given the GMM
ndim = size(data, 1);
C = sum(mu.*mu./sigma) + sum(log(sigma));
D = (1./sigma)' * (data .* data) - 2 * (mu./sigma)' * data  + ndim * log(2 * pi);
logprob = -0.5 * (bsxfun(@plus, C',  D));
logprob = bsxfun(@plus, logprob, log(w));

function y = logsumexp(x, dim)
% compute log(sum(exp(x),dim)) while avoiding numerical underflow
xmax = max(x, [], dim);
y    = xmax + log(sum(exp(bsxfun(@minus, x, xmax)), dim));
ind  = find(~isfinite(xmax));
if ~isempty(ind)
    y(ind) = xmax(ind);
end

F = reshape(F, ndim * nmix, 1); % F的维度:(ndim,nmix)
m = reshape(ubm.mu, ndim * nmix, 1); % 把均值列起来
F = F - N(idx_sv) .* m; % 消除一次充分统计量

EM算法迭代

EM目标优化公式：

$\begin{array}{l}{Q\left(T | T^{(t)}\right)=\sum_{s=1}^{S} E\left(\widehat{w}^{(t)}(s)^{T} T^{T} \Sigma^{-1} F(s)\right)} {-\frac{1}{2} \sum_{s=1}^{S} E\left(\widehat{w}^{(t)}(s)^{T} T^{T} \Sigma^{-1} N(s) T \widehat{w}^{(t)}(s)^{T}\right)}\end{array}$

E-Step：

EM算法的本质就是对隐变量进行估计，核心是隐变量的后验概率公式：

$y(s)|x \sim N(l^{-1}T^T\Sigma^{-1}F(s),l^{-1}(s))$

其中方差$l^{-1}(s)$为：

$\sigma(s) = l^{-1}(s) = I + T^T\Sigma^{-1}N(s)T$

有了上面的后验概率分布，我们可以确定新的"B-W"公式：

$\begin{aligned} E\left(\widehat{w}^{(t)}(s)^{T}\right)= & \sigma^{t}(s)^{-1} T^{T} \Sigma^{-1} F(s) \\ E\left(\widehat{w}^{(t)}(s)^{T} \cdot \widehat{w}^{(t)}(s)\right)= & \sigma^{t}(s)^{-1}+E\left(\widehat{w}^{(t)}(s)^{T}\right) \cdot E\left(\widehat{w}^{(t)}(s)^{T}\right)^{T} \end{aligned}$

说明：上面的第一个公式就是后验概率的均值

说明2： 第二个方差如何计算？ $E[X\cdot Y] = cov(X,Y)+E[X]E[Y]$

M-Step：

$T^{t+1}=argmaxQ(T|T^{(t)})$

因此有公式：

$\sum_{s=1}^S \Sigma^{-1}F(s)E(w^{(t)}(s)^T) = \sum_{s=1}^S \Sigma^{-1}N(s)T^{t+1}E(w^{(t)}(s)^T \cdot w^{(t)}(s))$

但是上面的方程没办法编程的原因在于$T$矩阵被夹在两个矩阵的中间，根本无法求逆，所以必须根据在分布内进行更新

上面的公式同上有：

$\sum_{s=1}^S F_c(s)E(w^{(t)}(s)^T) = T^{t+1}\sum_{s=1}^SN_c(s)E(w^{(t)}(s)^T \cdot w^{(t)}(s)) \\ T^{t+1} = (\sum_{s=1}^S F_c(s)E(w^{(t)}(s)^T))(\sum_{s=1}^SN_c(s)E(w^{(t)}(s)^T \cdot w^{(t)}(s)))^{-1}$

上面的方程与$Matlab$代码中的稍有不同，在于$T$矩阵的大小不同。

算法流程

1. 对$T$ 进行随机初始化

2. 根据$T^{(t)}$计算 $E\left(\widehat{w}^{(t)}(s)^{T}\right)$ 和 $E\left(\widehat{w}^{(t)}(s)^{T} \cdot \widehat{w}^{(t)}(s)\right)$

3. 重新计算$T^{(t+1)}$

4. 重复迭代第2,3步，大约5到10步迭代结束

5. 计算$i-vector$ ,即$E\left(\widehat{w}^{(t)}(s)^{T}\right)$

   $w = \sigma^{t}(s)^{-1} T^{T} \Sigma^{-1} F(s)$

Matlab代码编程

数据准备： stats(:)

stats = cell(nTra_Speakers, nTra_Channels);
for s=1:nTra_Speakers
    for c=1:nTra_Channels
        [N,F] = compute_bw_stats(trainSpeakerData{s,c}, ubm);
        stats{s,c} = [N; F]; %448 x 1？=[(32,1);(13*32,1)]
    end
end

充分统计量是两个列向量的叠加，N是$(nmix,1)$，F是超向量：$(nmix*ndim)$

核心主函数：train_tv_space

tvDim = 100;
niter = 5; % 迭代次数一般在5~10次。

%T矩阵(100,416）,之后还需要与i-vector相乘缩小维度
T = train_tv_space(stats(:), ubm, tvDim, niter, nworkers);

Function说明：

function T = train_tv_space(dataList, ubmFilename, tv_dim, niter, nworkers, tvFilename)
[ndim, nmix] = size(ubm.mu);
S = reshape(ubm.sigma, ndim * nmix, 1);

[N, F] = load_data(dataList, ndim, nmix); %dataList = stats(:)

% 初始化T矩阵
T = randn(tv_dim, ndim * nmix) * sum(S) * 0.001; 

fprintf('Re-estimating the total subspace with %d factors ...\n', tv_dim);

for iter = 1 : niter,
    fprintf('EM iter#: %d \t', iter);
    tim = tic;
    % E-step
    [LU, RU] = expectation_tv(T, N, F, S, tv_dim, nmix, ndim, nworkers);
    % M-step
    T = maximization_tv(LU, RU, ndim, nmix);
    tim = toc(tim);
    fprintf('[elaps = %.2f s]\n', tim);
end

1. 需要对$dataList$进行分析，划分相应的$N$和$F$

   stats{nTra_Speakers,channel}

   观察stats(:)，Matlab默认竖直读取数据。

   nfiles = length(datalist); % nfiles是总的语音文件 N = zeros(nfiles, nmix, 'single'); F = zeros(nfiles, ndim * nmix, 'single'); for file = 1 : nfiles, % datalist{file} 是stats文件，前1:nmix是N ，前nmix+1：end是F N(file, :) = datalist{file}(1:nmix); F(file, :) = datalist{file}(nmix + 1 : end); end

2. E - Step⭐️

   • 公式1：

   $\sigma(s) = l^{-1}(s) = I + T^T\Sigma^{-1}N(s)T$

   其中$\Sigma$ 是$KD\times KD$的块对角矩阵，子对角块是$\Sigma_1,\dots,\Sigma_C$ ;$N(s)$也是一个$KD\times KD$的块对角矩阵，子对角块是$N_1I,\dots,N_CI$。

   维度说明：$\sigma(s)$是隐变量的方差，故维度是

   $(tv,tv)=(tv,tv)+(KD,tv)'(KD,KD)*(KD,KD)*(KD,tv)$

   Matlab编程的方程为：

$\sigma(s) = I + (T./\Sigma).*N*T'$

其中：$\Sigma =\left[ \begin{matrix} \Sigma _1& & \\ & \ddots& \\ & & \Sigma _C\\ \end{matrix} \right]$ 简化为列矩阵，即$\Sigma =\left[ \begin{array}{c} \Sigma _1\\ \vdots\\ \Sigma _C\\ \end{array} \right]$

其中：$N\ =\ \left[ \begin{matrix} N_1I& & \\ & \ddots& \\ & & N_CI\\ \end{matrix} \right]$简化为列矩阵，即$N_2=\left[ \begin{array}{c} N_1I\\ \vdots\\ N_CI\\ \end{array} \right]$ ，为了进一步节省内存，可以简化为$N_1\ =\ \left[ \begin{array}{c} N_1\\ \vdots\\ N_C\\ \end{array} \right]$

idx_sv = reshape(repmat(1 : nmix, ndim, 1), ndim * nmix, 1)
   % 编程技巧：Matlab编程存储矩阵为N1，实际方程汇总为N2，故这里通过一串索引来扩展矩阵 N1 --> N2 
   N2 = N1(ix, idx_sv)

对于编程矩阵的简化：

$\begin{aligned} T\_invS = T^T\Sigma^{-1}=& \left[ \begin{matrix} a_{11}& \cdots& a_{1,KD}\\ \vdots& \ddots& \vdots\\ a_{tv\times 1}& \cdots& a_{tv\times \text{KD}}\\ \end{matrix} \right] _{tv,KD}\cdot \left[ \begin{matrix} \Sigma _1& & \\ & \ddots& \\ & & \Sigma _C\\ \end{matrix} \right] _{KD,KD} \\ =&\left[ \begin{matrix} \Sigma _1\cdot a_{11}& \cdots& \Sigma _1\cdot a_{1,KD}\\ \vdots& \ddots& \vdots\\ \Sigma _C\cdot a_{tv\times 1}& \cdots& \Sigma _C\cdot a_{tv\text{timesKD}}\\ \end{matrix} \right] \ \\ = &\ repmat\left( \left[ \begin{array}{c} \Sigma _1\\ \vdots\\ \Sigma _C\\ \end{array} \right] ,\left( 1,KD \right) \right) .*T \\ =& \ T./\Sigma 【维度为(tv,KD)】 \end{aligned}$

可以发现$T$矩阵乘以对角矩阵，相当于给T矩阵的每一个行向量乘上一个系数。

T_invS =  bsxfun(@rdivide, T, S');

同理也有下面：

$T\_invS_{(tv,KD)} * \left[ \begin{matrix} N_1I& & \\ & \ddots& \\ & & N_CI\\ \end{matrix} \right]_{KD,KD} =T\_invS.*repmat(\left[ \begin{array}{c} \N_1I\\ \vdots\\ \N _CI\\ \end{array} \right]_{KD,1},(1,tv))$

bsxfun(@times, T_invS, N1(ix, idx_sv)

综上，公式1的的Matlab编程：

$\sigma(s) = I + (T./\Sigma).*N*T'$

I = eye(tv_dim);
   L = I +  bsxfun(@times, T_invS, N1(ix, idx_sv)) * T';

期望公式：

$\begin{aligned} E\left(\widehat{w}^{(t)}(s)^{T}\right)= & \sigma^{t}(s)^{-1} T^{T} \Sigma^{-1} F(s) \\ E\left(\widehat{w}^{(t)}(s)^{T} \cdot \widehat{w}^{(t)}(s)\right)= & \sigma^{t}(s)^{-1}+E\left(\widehat{w}^{(t)}(s)^{T}\right) \cdot E\left(\widehat{w}^{(t)}(s)^{T}\right)^{T} \end{aligned}$

这里为了对每个高斯分布进行统计量操作，所以需要遍历所有的样本个数，将所有的Ex和Exx累加起来。

未批处理操作：

for (ix = 1 : size(F, 1), nworkers) % 这个是后验概率 y|x 的协方差的逆 L = I + bsxfun(@times, T_invS, N1(ix, idx_sv)) * T'; % this is the posterior covariance Cov(x,x) % 这个就是后验概率的协方差 Cxx = pinv(L); B = T_invS * F1(ix, :)'; Ex(:, ix) = Cxx * B; % 后验均值 E[w] (tv,1) Exx(:, :, ix) = Cxx + Ex(:, ix) * Ex(:, ix)';% 后验方差E[w^T*w],(tv,tv) end RU = RU + Ex * F1; % 因为这是批操作，所以要累加 for (mix = 1 : nmix, nworkers) tmp = bsxfun(@times, Exx, reshape(N1(:, mix),[1 1 len])); LU{mix} = LU{mix} + sum(tmp, 3); end

⭐️Matlab累加操作，是通过矩阵乘的方式完成的！！！！！

$Ex(w): (tv,nfiles)$

$Exx(w^Tw):(tv,tv,nfiles)$

$F:(KD,nfiles)$

⭐️累加操作：$\sum_{s=1}^S$

$Ex(w)*F'$以及$sum(Exx,3)$

完整代码如下：

function [LU, RU] = expectation_tv(T, N, F, S, tv_dim, nmix, ndim, nworkers)
   % compute the posterior means and covariance matrices of the factors 
   % or latent variables
   % T:(tv_dim,ndim * nmix)
   % N:(nfiles,nmix)
   % F:(nfiles,ndim * nmix)
   % S:(ndim * nmix, 1) -- ubm.Sigma
   
   idx_sv = reshape(repmat(1 : nmix, ndim, 1), ndim * nmix, 1);
   nfiles = size(F, 1);
   
   LU = cell(nmix, 1);
   LU(:) = {zeros(tv_dim)};
   RU = zeros(tv_dim, nmix * ndim);
   
   I = eye(tv_dim);
   T_invS =  bsxfun(@rdivide, T, S');
   
   parts = 250; % modify this based on your resources
   nbatch = floor( nfiles/parts + 0.99999 );
   for batch = 1 : nbatch,
       start = 1 + ( batch - 1 ) * parts;
       fin = min(batch * parts, nfiles);
       len = fin - start + 1;
       index = start : fin;
       N1 = N(index, :);
       F1 = F(index, :);
       Ex = zeros(tv_dim, len);
       Exx = zeros(tv_dim, tv_dim, len);
       parfor (ix = 1 : len, nworkers)
           % 这个是后验概率 y|x 的协方差的逆
           L = I +  bsxfun(@times, T_invS, N1(ix, idx_sv)) * T';  
           % this is the posterior covariance Cov(x,x) 
           % 这个就是后验概率的协方差
           Cxx = pinv(L);
           B = T_invS * F1(ix, :)';
           Ex(:, ix) = Cxx * B; % 后验均值 E[w] (tv,1)
           Exx(:, :, ix) = Cxx + Ex(:, ix) * Ex(:, ix)';% 后验方差E[w^T*w],(tv,tv)
       end
       RU = RU + Ex * F1; % 因为这是批操作，所以要累加
       parfor (mix = 1 : nmix, nworkers)
           tmp = bsxfun(@times, Exx, reshape(N1(:, mix),[1 1 len]));
           LU{mix} = LU{mix} + sum(tmp, 3);
       end
   end

3. M - Step

   $T$矩阵的更新公式：

   $T^{t+1} = (\sum_{s=1}^S F_c(s)E(w^{(t)}(s)^T))(\sum_{s=1}^SN_c(s)E(w^{(t)}(s)^T$

   - 代码中$RU$为$\sum_{s=1}^S F(s)E(w^{(t)}(s)^T)$ 可以直接计算出，大小为$(tv,KD)$，上面的等式中F带C，更新$T$矩阵的时候，只需要取出当前高斯分布下的D维，大小为$(tv,D)$ - 代码中$LU$ 为 $\sum_{s=1}^SN_c(s)E(w^{(t)}(s)^T \cdot w^{(t)}(s)))^{-1}$ 只有根据分布来计算，所以大小为$(nmix,1)$cell格式，每个数据的大小为$(tv,tv)$ - $T$ 矩阵的分块更新方式为：$T = LU$ \ $RU$ 维度的大小为【$(tv,D) = (tv,tv) *(tv,D)$】

   function RU = maximization_tv(LU, RU, ndim, nmix) % ML re-estimation of the total subspace matrix or the factor loading matrix for mix = 1 : nmix idx = ( mix - 1 ) * ndim + 1 : mix * ndim; %这里等式左边应该是T，但是为了节约内存所以直接调用RU这个变量了 RU(:, idx) = LU{mix} \ RU(:, idx); end

以上内容来源：https://blog.csdn.net/JackyTintin/article/details/79803501

模型理解：概率数据生成背景

PLDA是LDA的概率版本，两者都是将数据之间的差异分为类内差异和类间差异，但是从概率的角度进行建模。

定义：隐变量$y$决定每个高斯分布的生成方式，其均值的差异可以看做是类间方差$\Phi_b$ 。每个高斯分布的内部的数据差异可看做是类内差异$\Phi_w$。因此有以下的概率关系：

$\begin{array}{l} {p(\mathbf{x} | \mathbf{y})=\mathcal{N}\left(\mathbf{x} | \mathbf{y}, \Phi_{w}\right)} \\ {p(\mathbf{y})=\mathcal{N}\left(\mathbf{y} | \mathbf{m}, \Phi_{b}\right)} \end{array}$

LDA 假设各类中心服从离散分布，离散中心的个数固定，是训练数据中已知的类别数；PLDA 假设各类中心服从一个连续分布（高斯分布）。因此，PLDA 能够扩展到未知类别，从而用于未知类别的识别与认证。

这里要求协方差矩阵 $\Phi_w$ 是正定的对称方阵，反映了类内（within-class）的差异；$\Phi_b$ 是半正定 的对称方阵，反映了类间（between-class）的差异。因此，PLDA 建模了数据生成的过程，并且同时 LDA 一样，显式地考虑了类内和类间方差。

为了推导方便，消除类内方差为单位矩阵：存在一个矩阵$V$可以让$\Phi_w$ 和$\Phi_b$同时合同对角化。存在可逆矩阵$V$,满足$V^T\Phi_bV =\Psi$ (对角阵) 且 $V^T\Phi_wV=I$（单位对角矩阵）。

PLDA等价表述：

$x = m + A u$

其中，

$u \sim N(\cdot|v,I) \\ v \sim N(\cdot |0,\Psi) \\ A = (V^T)^{-1}$

$u$ 是数据空间在投影空间的对应投影点。$\Psi$反映了类内（$within-class$）的差异； $I$ 反映了类间（$between-class$）的差异（这里被规一化为单位矩阵）。

基于PLDA的推断

对于每一个观测数据 $\mathbf{x}$ 我们都可以计算对应的 $u=A^{-1}(x-m)$。$PLDA$的推断的投影空间中进行

给定观一组同类的测数据 $\mathbf{u}_{1,\dots,n} ,\mathbf{v}$的后验概率分布为：

$p\left(\mathbf{v} | \mathbf{u}_{1, \ldots, n}\right)=\mathcal{N}\left(\mathbf{v} | \frac{n \Psi}{n \Psi+I} \overline{\mathbf{u}}, \frac{\Psi}{n \Psi+I}\right)$

如何求解隐变量的后验概率？

给定某类的 $n$ 个数据 $x_{1,\dots,n}$,则$y$的后验分布可以为：

$p\left(y | x_{1, \ldots, n}\right)=p\left(x_{1, \ldots, n} | y\right) p(y) / p\left(x_{1, \ldots, n}\right) \propto p\left(x_{1, \ldots, n} | y\right) p(y)$

后验为两个高斯分布的乘积，因此也服从高斯。因此，我们只需要计算均值向量和方差矩阵，即可以确定后验分布。

$\begin{aligned} \ln p\left(y | x_{1, \dots, n}\right)= &\ln p\left(x_{1}, \ldots, n | y\right)+\ln p(y) \\ = &\sum_{i} \ln p\left(x_{i} | y\right)+\ln p(y)\\ = &C-0.5 \sum_{i} y^{T} \Phi_{w}^{-1} y+\sum_{i} x_{i}^{T} \Phi_{w}^{-1} y-0.5 y^{T} \Phi_{b}^{-1} y+m^{T} \Phi_{b}^{-1} y \\ = & C - 0.5y^T(n\Phi_w^{-1}+\Phi_b^{-1})y+\sum_{i} x_{i}^{T} \Phi_{w}^{-1}y+m^T\Phi_b^{-1}y \end{aligned}$

以上第三个等式的来源：

$p(x|y) \sim N(x|y,\Phi_w) \\ p(y) \sim N(y|m,\Phi_b)$

高斯分布公式为：

$\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \mathbf{\Sigma})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{T} \mathbf{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\}$

展开后，忽略常数项，只关注指数项里面的系数：

$\begin{aligned} & exp\left\{ -0.5x^T\Sigma^{-1}x + x^T \Sigma^{-1}\mu \right\} \\ =&exp\left\{ -0.5x^T\Sigma^{-1}x + \mu ^T(\Sigma^{-1})^T x \right\} \\ =&exp\left\{ -0.5x^T\Sigma^{-1}x + \mu ^T\Sigma^{-1} x \right\} \end{aligned}$

根据高斯分布的一次项和 二次项可得 均值 和 方差 ：

$\begin{aligned} \hat{m} =& \hat{\Phi}(\Phi_b^{-1}m+\Phi_w^{-1}\sum_ix_i) \\ \hat{\Phi} =& (n\Phi_w^{-1}+\Phi_b^{-1})^{-1} \\ \end{aligned}$

综上后验分布为：

$p(y|x_1,\dots,n) \sim N(\hat{m},\hat{\Phi})$

PLDA 在声纹识别中的判断

因此，对于未知数据点$u^p$ 以及某类的若干数据点 $\mathbf{u}^g_{1,\dots,n}$ 属于该类的似然值：

用同样的方法，可以得到下面的式子：

$\begin{array}{c} {p\left(\mathbf{u}^{p} | \mathbf{u}_{1, \ldots, n}^{g}\right)=\mathcal{N}\left(\frac{n \Psi}{n \Psi+I} \overline{\mathbf{u}}^{g}, \frac{\Psi}{n \Psi+I}+I\right)} \\ {\ln p\left(\mathbf{u}^{p} | \mathbf{u}_{1, \ldots, n}^{g}\right)=C-\frac{1}{2}\left(\mathbf{u}^{p}-\frac{n \Psi}{n \Psi+I} \overline{\mathbf{u}}^{g}\right)^{T}\left(\frac{\Psi}{n \Psi+I}+I\right)^{-1}\left(\mathbf{u}^{p}-\frac{n \Psi}{n \Psi+I} \overline{\mathbf{u}}^{g}\right)-\frac{1}{2} \ln \left|\frac{\Psi}{n \Psi+I}+I\right|} \end{array}$

其中，$C = -\frac{1}{2}dln2\pi$ 是与数据无关的常量，$d$是 数据的维度。特殊的，$u^p$ 不属于任何已知类的概率为：

$\begin{aligned} & p\left(\mathbf{u}^{p} | \phi\right)=\mathcal{N}\left(\mathbf{u}^{p} | 0, \Psi+I\right) \\ \ln p\left(\mathbf{u}^{p} | 0\right)=& C-\frac{1}{2} \mathbf{u}^{p T}(\Psi+I)^{-1} \mathbf{u}^{p}-\frac{1}{2} \ln |\Psi+I| \end{aligned}$

由于 $\Psi$ 是对角阵，因此上式中各个协方差也都是对角阵，因此，似然和对数似然都很容易求得。

利用 PLDA 进行识别（recognition）方法如下：

$i=\underset{i}{\operatorname{argmax}} \ln p\left(\mathbf{u}^{p} | \mathbf{u}_{1, \ldots, n}^{g_{i}}\right)$

对于认证问题（verification），可以计算其似然比：

$R = \frac{p(u^p|u^g_{1,\dots,n})}{p(u^p|\phi)}$

或似然比对数（log likelihood ratio）：

$lnR = lnp(u^p|u^g_{1,\dots,n})-lnp(u^p|\phi)$

适当的选定阈值 $T$，当 $R>T$ 判定 $\mathbf{u}$ 与已知数据属于同一个类别，反之则不是。

期望最大化方法（EM）

PLDA 中，我们需要估计的参数包括 $A、\Psi 和 \mathbf{m}$。

算法流程

输入：$K$ 类 $d$ 维数据，第 $k$ 个类别包含 $n_k$ 个样本，记 $x_{ki}∈R^d,1≤k≤K$为 第 $k$ 个类别的第 $i$ 个样本点。

输出：$d\times d$对称矩阵$\Phi_w$，$d \times d$ 对称矩阵$\Phi_b$,$d$ 维向量 $m$。

1. 计算统计量

   $\begin{aligned} N=&\sum_{k=1}^{K} n_{k}, \\ f_{k}=&\sum_{i=1}^{n_{k}} x_{k i}\\ m=&\frac{1}{N} \sum_{k} f_{k}\\ S=&\sum_{k} \sum_{i} x_{k i} x_{k i}^{T} \end{aligned}$

2. 随机初始化$\Phi_w,\Phi_b,m$

3. 重复如下步骤至到满足终止条件：

   1. 对每一个类别，计算：$\hat{\Phi}=\left(n \Phi_{w}^{-1}+\Phi_{b}^{-1}\right)^{-1}, y=\hat{\Phi}\left(\Phi_{b}^{-1} m+\Phi_{w}^{-1} f\right), y y t=\hat{\Phi}+y y^{T}$
   2. 中间计算：$R=\sum_{k} n_{k} \cdot y y t_{k}, T=\sum_{k} y_{k} f_{k}^{T}, P=\sum_{k} \hat{\Phi}_{k}, \quad E=\sum_{k}\left(y_{k}-m\right)\left(y_{k}-m\right)^{T}$
   3. 更新：$\Phi_{w}=\frac{1}{N}\left(S+R-\left(T+T^{T}\right)\right), \Phi_{b}=\frac{1}{K}(P+E)$

基于$\Phi_w $和$\Phi_b$可以计算出$\Psi$和$A^{-1}$。

根据LDA的结论，最大投影方向为 ：$\Phi_w^{-1}\Phi_b$的特征向量$w_1,\dots,d$,每个特征向量为一列，组成矩阵 $W$。则有：

$\begin{aligned} &\Lambda_{b}=W^{T} \Phi_{b} W\\ &\Lambda_{w}=W^{T} \Phi_{w} W \end{aligned}$

为了让$\Phi_w$合同对角化到单位矩阵，显然有下式：

$I=\Lambda_{w}^{-1 / 2} \Lambda_{w} \Lambda_{w}^{-1 / 2}=\Lambda_{w}^{-1 / 2} W^{T} \Phi_{w} W \Lambda_{w}^{-1 / 2}$

则找到了转换矩阵$V$：

$V = W \Lambda_{w}^{-1 / 2}$

$\begin{aligned} \Psi =& V^T\Phi_bV=(W \Lambda_{w}^{-1 / 2})^T\Phi_bW \Lambda_{w}^{-1 / 2}\\ = &(\Lambda_{w}^{-1 / 2})^T\Lambda_b \Lambda_{w}^{-1 / 2} \\ = &上面都是对称矩阵所以相当于点乘 \\ = &\Lambda_b \Lambda_w^{-1} \end{aligned}$

至于A的求解方法，还记得PLDA的等效公式吗？：

$x = m + A u$

$u$是属于投影空间，$x$属于实际空间，故说明$A$是连接实际空间与投影空间之间的变换矩阵，故：

$A = V^{-1}$

PLDA算法版本2 – MSR工具箱

数据生成背景：

$\mathbf{x}_{i j}=\mu+\mathbf{F h}_{i}+\mathbf{G} \mathbf{w}_{i j}+\epsilon_{i j}$

上述超参数说明：

在文献【1】中将 $F$ 和 $G$ 合并起来：

$\left[\begin{array}{c} {\mathbf{x}_{1}} \\ {\mathbf{x}_{2}} \\ {\vdots} \\ {\mathbf{x}_{N}} \end{array}\right]=\left[\begin{array}{c} {\mu} \\ {\mu} \\ {\vdots} \\ {\mu} \end{array}\right]+\left[\begin{array}{ccccc} {\mathbf{F}} & {\mathbf{G}} & {\mathbf{0}} & {\ldots} & {\mathbf{0}} \\ {\mathbf{F}} & {\mathbf{0}} & {\mathbf{G}} & {\ldots} & {\mathbf{0}} \\ {\vdots} & {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {\mathbf{F}} & {\mathbf{0}} & {\mathbf{0}} & {\ldots} & {\mathbf{G}} \end{array}\right]\left[\begin{array}{c} {\mathbf{h}} \\ {\mathbf{w}_{1}} \\ {\mathbf{w}_{2}} \\ {\vdots} \\ {\mathbf{\epsilon}_{N}} \end{array}\right]+\left[\begin{array}{c} {\epsilon_{1}} \\ {\epsilon_{2}} \\ {\vdots} \\ {\epsilon_{N}} \end{array}\right]$

简化为：

$\mathbf{x}^{\prime}=\mu^{\prime}+\mathbf{A} \mathbf{y}+\epsilon^{\prime}$

E步公式：

$\begin{aligned} E\left[\mathbf{y}_{i}\right] &=\left(\mathbf{A}^{T} \Sigma^{\prime-1} \mathbf{A}+\mathbf{I}\right)^{-1} \mathbf{A}^{T} \Sigma^{\prime-1}\left(\mathbf{x}_{i}-\mu^{\prime}\right) \\ E\left[\mathbf{y}_{i} \mathbf{y}_{i}^{T}\right] &=\left(\mathbf{A}^{T} \Sigma^{\prime-1} \mathbf{A}^{T}+\mathbf{I}\right)^{-1}+E\left[\mathbf{y}_{i}\right] E\left[\mathbf{y}_{i}\right]^{T} \end{aligned}$

M步公式：

$\begin{array}{l} {\mu=\frac{1}{I J} \sum_{i, j} \mathbf{x}_{i j}} \\ {\mathbf{B}=\left(\sum_{i, j}\left(\mathbf{x}_{i j}-\mu\right) E\left[\mathbf{z}_{i}\right]^{T}\right)\left(\sum_{i, j} E\left[\mathbf{z}_{i} \mathbf{z}_{i}^{T}\right]\right)^{-1}} \\ {\mathbf{\Sigma}=\frac{1}{I J} \sum_{i, j} \operatorname{Diag}\left[\left(\mathbf{x}_{i j}-\mu\right)\left(\mathbf{x}_{i j}-\mu\right)^{T}-\mathbf{B} E\left[\mathbf{z}_{i}\right]\left(\mathbf{x}_{i j}-\mu\right)^{T}\right]} \end{array}$

在文献【2】讲述了识别原理：

对于two-covariance版本，需要同时迭代出两个方差，文献【1】中采用合并的方式训练，上述模型中包含两个部分：

• (i) 【类间】说话人的部分：$\mathbf{x}_{i j}=\mu+\mathbf{F h}_{i}$
• (ii) 【类内】通道效应部分：$c_r = \mathbf{G} \mathbf{w}_{i j}+\epsilon_{i j}$

而这里我们只对说话人的进行建模，而对通道效应进行简化：

$\mathbf{x}_{i j} =\mu+\mathbf{F h}_{i}+\epsilon_{i j}$

这里$\epsilon_{i j}$从原先的diagonal covariace$\Sigma$ 变为full-covariance $\Sigma$

Verification score

算法原理

得分公式如下：

$\text { score }=\log \frac{p\left(\eta_{1}, \eta_{2} | \mathscr{H}_{s}\right)}{p\left(\eta_{1} | \mathscr{H}_{d}\right) p\left(\eta_{2} | \mathscr{H}_{d}\right)}$

分子代表两个不同的i-vector来源于同一个分布的概率；而分母代表i-vector来源于不同的分布的概率。分子应越大越好，分母应越小越好，得分越大说明两者为同一人的概率越高。

如果有PPCA或者因子分析法的推导基础，很容易可以写出两个变量的似然概率函数，来源于同一个分布，方差应为fuu-covariance而不同分布的变量独立，所以为diagnonal-covariance。

$\begin{aligned} \text { score } &=\log \mathcal{N}\left(\left[\begin{array}{l} \eta_{1} \\ \eta_{2} \end{array}\right] ;\left[\begin{array}{l} m \\ m \end{array}\right],\left[\begin{array}{ll} \Sigma_{t o t} & \Sigma_{a c} \\ \Sigma_{a c} & \Sigma_{t o c} \end{array}\right]\right) \\ &-\log \mathcal{N}\left(\left[\begin{array}{l} \eta_{1} \\ \eta_{2} \end{array}\right] ;\left[\begin{array}{l} m \\ m \end{array}\right],\left[\begin{array}{cc} \Sigma_{t o t} & 0 \\ 0 & \Sigma_{t o t} \end{array}\right]\right) \end{aligned}$

这里文献【2】中写错了，其中$\Sigma_{tot}=FF^T+\Sigma$ 和$\Sigma_{ac} = FF^T$。

文献【3】给出了有关矩阵操作的所有数学知识点，可以随时翻阅，以下是第9章的Special Matrix的内容

• 这里需要用到分块矩阵的协方差

$\begin{aligned}\left[\begin{array}{c|c}\mathbf{A}_{11} & \mathbf{A}_{12} \\\hline \mathbf{A}_{21} & \mathbf{A}_{22}\end{array}\right]^{-1}= &\left[\begin{array}{c|c}\mathbf{C}_{1}^{-1} & -\mathbf{A}_{11}^{-1} \mathbf{A}_{12} \mathbf{C}_{2}^{-1} \\\hline-\mathbf{C}_{2}^{-1} \mathbf{A}_{21} \mathbf{A}_{11}^{-1} & \mathbf{C}_{2}^{-1}\end{array}\right] \\= &\left[\begin{array}{c|c}\mathbf{A}_{11}^{-1}+\mathbf{A}_{11}^{-1} \mathbf{A}_{12} \mathbf{C}_{2}^{-1} \mathbf{A}_{21} \mathbf{A}_{11}^{-1} & -\mathbf{C}_{1}^{-1} \mathbf{A}_{12} \mathbf{A}_{22}^{-1} \\\hline-\mathbf{A}_{22}^{-1} \mathbf{A}_{21} \mathbf{C}_{1}^{-1} & \mathbf{A}_{22}^{-1}+\mathbf{A}_{22}^{-1} \mathbf{A}_{21} \mathbf{C}_{1}^{-1} \mathbf{A}_{12} \mathbf{A}_{22}^{-1}\end{array}\right]\end{aligned}$

• 其中：

$\begin{array}{l}\mathbf{C}_{1}=\mathbf{A}_{11}-\mathbf{A}_{12} \mathbf{A}_{22}^{-1} \mathbf{A}_{21} \\\mathbf{C}_{2}=\mathbf{A}_{22}-\mathbf{A}_{21} \mathbf{A}_{11}^{-1} \mathbf{A}_{12}\end{array}$

根据上述注释中的知识点，可以推导出：

$score = \eta_1^TQ\eta + \eta_2^TQ\eta_2+2\eta_1^TP\eta_1+const$

其中：

$\begin{array}{l} \mathrm{Q}=\Sigma_{\text {tot}}^{-1}-\left(\Sigma_{\text {tot}}-\Sigma_{a c} \Sigma_{\text {tot}}^{-1} \Sigma_{a c}\right)^{-1} \\ \mathrm{P}=\Sigma_{\text {tot}}^{-1} \Sigma_{a c}\left(\Sigma_{\text {tot}}-\Sigma_{a c} \Sigma_{\text {tot}}^{-1} \Sigma_{a c}\right)^{-1} \end{array}$

这里还可以进一步进行降维处理：

1. 将$P$矩阵进行特征对角化

$\begin{aligned} \mathrm{P} &=\left[\mathrm{U}_{K} | \mathrm{U}_{D-K}\right] \operatorname{diag}\left(\left[\lambda_{1}, \ldots, \lambda_{K}, 0,0, \ldots\right]\right)\left[\mathrm{U}_{K} | \mathrm{U}_{D-K}\right]^{T} \\ &=\mathrm{U}_{K} \operatorname{diag}\left(\left[\lambda_{1}, \ldots, \lambda_{K}\right]\right) U_{K}^{T} \end{aligned}$

2. 得到新的score公式，此时$\eta_{1}$从$(D,1)$维降维到$\tilde{\eta}_{1}$得$(K,1)$维

   $\text { score }=\tilde{\eta}_{1}^{T} \widetilde{\mathrm{Q}} \tilde{\eta}_{1}+\tilde{\eta}_{2}^{T} \widetilde{\mathrm{Q}} \tilde{\eta}_{2}+2 \tilde{\eta}_{1}^{T} \Lambda \tilde{\eta}_{2}+\text { const }$

   其中$\Lambda = diag([\lambda_1,\dots,\lambda_K])$。

   这样得做法的目的就为了使得投影后的$\eta_1$和$\eta_2$的协方差矩阵变成diagonal covariance，从而达到了消除了不同人之间得相关性（即类间方差）的目的。

Matlab中MSR中的得分内容

function scores = score_gplda_trials(plda, model_iv, test_iv)

Phi     = plda.Phi;
Sigma   = plda.Sigma;
W       = plda.W;
M       = plda.M;

%%%%% post-processing the model i-vectors %%%%%
model_iv = bsxfun(@minus, model_iv, M); % centering the data
model_iv = length_norm(model_iv); % normalizing the length
model_iv = W' * model_iv; % whitening data

%%%%% post-processing the test i-vectors %%%%%
test_iv = bsxfun(@minus, test_iv, M); % centering the data
test_iv = length_norm(test_iv); % normalizing the length
test_iv  = W' * test_iv; % whitening data

nphi = size(Phi, 1); 

% 根据论文的Section2.3 Verification score的公式
% 1. G-PLDA case的公式
Sigma_ac  = Phi * Phi';
Sigma_tot = Sigma_ac + Sigma;

Sigma_tot_i = pinv(Sigma_tot);
Sigma_i = pinv(Sigma_tot-Sigma_ac*Sigma_tot_i*Sigma_ac);
Q = Sigma_tot_i-Sigma_i;
P = (Sigma_tot_i * Sigma_ac) * Sigma_i;

% 2. 根据论文，可以进一步对P进行降维处理
[U, S] = svd(P);
S = diag(S);
Lambda = diag(S(1 : nphi));
Uk     = U(:, 1 : nphi);
Q_hat  = Uk' * Q * Uk;

model_iv = Uk' * model_iv;
test_iv  = Uk' * test_iv;

score_h1 = diag(model_iv' * Q_hat * model_iv);
score_h2 = diag(test_iv' * Q_hat * test_iv);
score_h1h2 = 2 * model_iv' * Lambda * test_iv;

scores = bsxfun(@plus, score_h1h2, score_h1);
scores = bsxfun(@plus, scores, score_h2');

文献：

1. Probabilistic Linear Discriminant Analysis for Inferences About Identity

2. Analysis of I-vector Length Normalization in Speaker Recognition Systems

3. The Matrix Cookbook

4. EM for Probabilistic LDA 详细推导了PLDA的EM算法（略复杂）

   https://sites.google.com/site/fastplda/ 提供了对应的代码