Kaldi中的ivector-原理和实现

# computing a PCA transform from the fbank data."
steps/online/nnet2/get_pca_transform.sh
# training the diagonal UBM."
steps/online/nnet2/train_diag_ubm.sh 
# training the iVector extractor"
steps/online/nnet2/train_ivector_extractor.sh

gmm_feats="ark,s,cs:apply-cmvn-online --config=$dir/online_cmvn.conf --spk2utt=ark:$sdata/JOB/spk2utt $dir/global_cmvn.stats scp:$sdata/JOB/feats.scp ark:- | splice-feats $splice_opts ark:- ark:- | transform-feats $dir/final.mat ark:- ark:- | subsample-feats --n=$subsample ark:- ark:- |"

feats="ark,s,cs:splice-feats $splice_opts scp:$sdata/JOB/feats.scp ark:- | transform-feats $dir/final.mat ark:- ark:- | subsample-feats --n=$subsample ark:- ark:- |"

ivector-extractor-init --ivector-dim=$ivector_dim --use-weights=false \
     "gmm-global-to-fgmm $dir/final.dubm -|" $dir/0.ie

gmm-global-get-post --n=$num_gselect --min-post=$min_post $dir/final.dubm "$gmm_feats" ark:- \| \
    scale-post ark:- $modified_posterior_scale "ark:|gzip -c >$dir/post.JOB.gz" || exit 1;

Args=() # bash array of training commands for 1:nj, that put accs to stdout.
for j in $(seq $nj_full); do
  Args[$j]=`echo "ivector-extractor-acc-stats --num-threads=$num_threads $dir/$x.ie '$feats' 'ark,s,cs:gunzip -c $dir/post.JOB.gz|' -|" | sed s/JOB/$j/g`
done

for g in $(seq $nj); do
  start=$[$num_processes*($g-1)+1]
  ivector-extractor-sum-accs --parallel=true "${Args[@]:$start:$num_processes}" $dir/acc.$x.$g 
done

for j in $(seq $nj); do
      accs+="$dir/acc.$x.$j "
done

ivector-extractor-sum-accs $accs $dir/acc.$x 
ivector-extractor-est --num-threads=$nt $dir/$x.ie $dir/acc.$x $dir/$[$x+1].ie 

ivector-extractor-acc-stats 
{
  utt_stats.AccStats(feats, post); //利用r计算0阶/1阶。
  CommitStatsForUtterance(extractor, utt_stats);//计算C/A，mean,Var
}

void IvectorExtractorConv::GetStats(
    const MatrixBase<BaseFloat> &feats,
    const MatrixBase<BaseFloat> &post,
    IvectorExtractorConvUtteranceStats *stats) const {
  typedef std::vector<std::pair<int32, BaseFloat> > VecType;

  int32 num_frames = feats.NumRows(), num_gauss = NumGauss(),
      feat_dim = FeatDim();
  KALDI_ASSERT(feats.NumCols() == feat_dim);
  KALDI_ASSERT(stats->gamma.Dim() == num_gauss &&
               stats->X.NumCols() == feat_dim);
  bool update_variance = (!stats->S.empty());
  
  for (int32 t = 0; t < num_frames; t++) {
    SubVector<BaseFloat> frame(feats, t);
    SpMatrix<double> outer_prod;
    if (update_variance) {
      outer_prod.Resize(feat_dim);
      outer_prod.AddVec2(1.0, frame);
    }
    for (int32 i = 0; i < num_gauss; i++) {
      double weight = post(t,i);
      stats->gamma(i) += weight;
      stats->X.Row(i).AddVec(weight, frame);
      if (update_variance)
        stats->S[i].AddSp(weight, outer_prod);
    }
  }

  for (int32 i = 0; i < num_gauss; i++) {
      stats->X.Row(i).AddVec(-stats->gamma(i), Means_.Row(i));
  }
}

CommitStatsForUtterance
{
 extractor.GetIvectorDistribution(utt_stats,
                                   &ivec_mean,
                                   &ivec_var);//mean,var
  CommitStatsForM(extractor, utt_stats, ivec_mean, ivec_var);//计算C/A，
}

GetIvectorDistribution
{
  GetIvectorDistMean(utt_stats, &linear, &quadratic);
    GetIvectorDistPrior(utt_stats, &linear, &quadratic);
    if (var != NULL) {
      var->CopyFromSp(quadratic);
      var->Invert(); // now it's a variance.

      // mean of distribution = quadratic^{-1} * linear...
      mean->AddSpVec(1.0, *var, linear, 0.0);
    } else {
      quadratic.Invert();
      mean->AddSpVec(1.0, quadratic, linear, 0.0);
    }
}

void IvectorExtractor::GetIvectorDistMean(
    const IvectorExtractorUtteranceStats &utt_stats,
    VectorBase<double> *linear,
    SpMatrix<double> *quadratic) const {
  int32 I = NumGauss();
  for (int32 i = 0; i < I; i++) {
    double gamma = utt_stats.gamma_(i);
    if (gamma != 0.0) {
      SubVector<double> x(utt_stats.X_, i); // == \gamma(i) \m_i
      // next line: a += \gamma_i \M_i^T \Sigma_i^{-1} \m_i
      linear->AddMatVec(1.0, Sigma_inv_M_[i], kTrans, x, 1.0);
    }
  }
  SubVector<double> q_vec(quadratic->Data(), IvectorDim()*(IvectorDim()+1)/2);
  q_vec.AddMatVec(1.0, U_, kTrans, utt_stats.gamma_, 1.0);
}

double IvectorExtractorStats::Update(
    const IvectorExtractorEstimationOptions &opts,
    IvectorExtractor *extractor) const {
  CheckDims(*extractor);
  if (tot_auxf_ != 0.0) {
    KALDI_LOG << "Overall auxf/frame on training data was "
              << (tot_auxf_/gamma_.Sum()) << " per frame over "
              << gamma_.Sum() << " frames.";
  }

  double ans = 0.0;
  ans += UpdateProjections(opts, extractor);
  if (extractor->IvectorDependentWeights())
    ans += UpdateWeights(opts, extractor);
  if (!S_.empty())
    ans += UpdateVariances(opts, extractor);
  ans += UpdatePrior(opts, extractor); // This will also transform the ivector
                                       // space.  Note: this must be done as the
                                       // last stage, because it will make the
                                       // stats invalid for that model.
  KALDI_LOG << "Overall objective-function improvement per frame was " << ans;
  extractor->ComputeDerivedVars();
  return ans;
}

// Note, throughout this file we use SGMM-type notation because
// that's what I'm comfortable with.
// Dimensions:
//  D is the feature dim (e.g. D = 60)
//  I is the number of Gaussians (e.g. I = 2048)
//  S is the ivector dim (e.g. S = 400)

  /// R_i, quadratic term for ivector subspace (M matrix)estimation.  This is a
  /// kind of scatter of ivectors of training speakers, weighted by count for
  /// each Gaussian.  Conceptually vector<SpMatrix<double> >, but we store each
  /// SpMatrix as a row of R_.  Conceptually, the dim is [I][S][S]; the actual
  /// dim is [I][S*(S+1)/2].
  Matrix<double> R_;

  /// Weight projection vectors, if used.  Dimension is [I][S]
  Matrix<double> w_;

  /// If we are not using weight-projection vectors, stores the Gaussian mixture
  /// weights from the UBM.  This does not affect the iVector; it is only useful
  /// as a way of making sure the log-probs are comparable between systems with
  /// and without weight projection matrices.
  Vector<double> w_vec_;

  /// Ivector-subspace projection matrices, dimension is [I][D][S].
  /// The I'th matrix projects from ivector-space to Gaussian mean.
  /// There is no mean offset to add-- we deal with it by having
  /// a prior with a nonzero mean.
  std::vector<Matrix<double> > M_;

void IvectorExtractorConv::GetIvectorDistribution(
    const IvectorExtractorConvUtteranceStats &utt_stats,
    VectorBase<double> *mean,
    SpMatrix<double> *var) const {
    Vector<double> linear(IvectorDim());
    SpMatrix<double> quadratic(IvectorDim());
    GetIvectorDistMean(utt_stats, &linear, &quadratic);
    if (var != NULL) {
      var->CopyFromSp(quadratic);
      var->Invert(); // now it's a variance.

      // mean of distribution = quadratic^{-1} * linear...
      mean->AddSpVec(1.0, *var, linear, 0.0);
    } else {
      quadratic.Invert();
      mean->AddSpVec(1.0, quadratic, linear, 0.0);
    }
}

void IvectorConvStats::CommitStatsForM(
    const IvectorExtractorConv &extractor,
    const IvectorExtractorConvUtteranceStats &utt_stats,
    const VectorBase<double> &ivec_mean,
    const SpMatrix<double> &ivec_var) {
  subspace_stats_lock_.lock();

  // We do the occupation stats here also.
  gamma_.AddVec(1.0, utt_stats.gamma);
  
  // Stats for the linear term in M:
  for  (int32 i = 0; i < extractor.NumGauss(); i++) {
    Y_[i].AddVecVec(1.0, utt_stats.X.Row(i),
                    Vector<double>(ivec_mean));
  }

  int32 ivector_dim = extractor.IvectorDim();
  // Stats for the quadratic term in M:
  SpMatrix<double> ivec_scatter(ivec_var);
  ivec_scatter.AddVec2(1.0, ivec_mean);
  SubVector<double> ivec_scatter_vec(ivec_scatter.Data(),
                                     ivector_dim * (ivector_dim + 1) / 2);
  R_.AddVecVec(1.0, utt_stats.gamma, ivec_scatter_vec);

  subspace_stats_lock_.unlock();
}

ivec_mean -> w
M->T
Y->C
R->A
U->T^T Sigma T
  /// U_i = M_i^T \Sigma_i^{-1} M_i
X_  -> f
i -> c

double impr = SolveQuadraticMatrixProblem(R, Y_[i], SigmaInv, solver_opts, &M),

double IvectorConvStats::UpdateProjection(
    const IvectorExtractorConvEstimationOptions &opts,
    int32 i,
    IvectorExtractorConv *extractor) const {
  int32 I = extractor->NumGauss(), S = extractor->IvectorDim();
  KALDI_ASSERT(i >= 0 && i < I);
  /*
    For Gaussian index i, maximize the auxiliary function
       Q_i(x) = tr(M_i^T Sigma_i^{-1} Y_i)  - 0.5 tr(Sigma_i^{-1} M_i R_i M_i^T)
   */
  if (gamma_(i) < opts.gaussian_min_count) {
    KALDI_WARN << "Skipping Gaussian index " << i << " because count "
               << gamma_(i) << " is below min-count.";
    return 0.0;
  }
  SpMatrix<double> R(S, kUndefined), SigmaInv(extractor->Sigma_inv_[i]);
  SubVector<double> R_vec(R_, i); // i'th row of R; vectorized form of SpMatrix.
  SubVector<double> R_sp(R.Data(), S * (S+1) / 2);
  R_sp.CopyFromVec(R_vec); // copy to SpMatrix's memory.

  Matrix<double> M(extractor->M_[i]);
  SolverOptions solver_opts;
  solver_opts.name = "M";
  solver_opts.diagonal_precondition = true;
  // TODO: check if inversion is sufficient?
  double impr = SolveQuadraticMatrixProblem(R, Y_[i], SigmaInv, solver_opts, &M),
      gamma = gamma_(i);
  if (i < 4) {
    KALDI_VLOG(1) << "Objf impr for M for Gaussian index " << i << " is "
                  << (impr / gamma) << " per frame over " << gamma << " frames.";
  }
  extractor->M_[i].CopyFromMat(M);
  return impr;
}

local/chain/train_ivector.sh: computing a PCA transform from the fbank data.
feat-to-dim scp:exp/chain/diag_ubm_fbank/train_subset/feats.scp -
steps/online/nnet2/get_pca_transform.sh --cmd queue.pl --splice-opts --left-context=3 --right-context=3 --max-utts 10000 --subsample 2 --dim 23 exp/chain/diag_ubm_fbank/train_subset exp/chain/pca_transform_fbank
Done estimating PCA transform in exp/chain/pca_transform_fbank

// bin/est-pca.cc
est-pca --dim=$dim --normalize-variance=$normalize_variance \
    --normalize-mean=$normalize_mean "$feats" $dir/final.mat || exit 1;

kaldi/src/bin/matrix-dim exp/chain/pca_transform_fbank/final.mat

/export/maryland/chaoyang/kaldi/src/bin/copy-matrix --binary=false exp/chain/pca_transform_fbank/final.mat -

steps/online/nnet2/train_diag_ubm.sh --cmd queue.pl --nj 15 --num-frames 700000 --num-threads 8 exp/chain/diag_ubm_fbank/train_subset 512 exp/chain/pca_transform_fbank exp/chain/diag_ubm_fbank

gmm-global-init-from-feats --num-threads=$num_threads --num-frames=$num_frames --min-gaussian-weight=$min_gaussian_weight \
    --num-gauss=$num_gauss --num-gauss-init=$num_gauss_init --num-iters=$num_iters_init \
    "$all_feats" $dir/0.dubm 

gmm-gselect --n=$num_gselect $dir/0.dubm "$feats" "ark:|gzip -c >$dir/gselect.JOB.gz";

for x in `seq 0 $[$num_iters-1]`; do
      gmm-global-acc-stats "--gselect=ark,s,cs:gunzip -c $dir/gselect.JOB.gz|" $dir/$x.dubm "$feats" $dir/$x.JOB.acc
      gmm-global-est $opt --min-gaussian-weight=$min_gaussian_weight $dir/$x.dubm "gmm-global-sum-accs - $dir/$x.*.acc|" $dir/$[$x+1].dubm
done
mv $dir/$num_iters.dubm $dir/final.dubm;

<DiagGMM>
<GCONSTS>[512]
<WEIGHTS>[512]
<MEANS_INVVARS>  [512 * 23]
<INV_VARS>[512 * 23]
</DiagGMM>

  /// Equals log(weight) - 0.5 * (log det(var) + mean*mean*inv(var))
  Vector<BaseFloat> gconsts_;
  bool valid_gconsts_;   ///< Recompute gconsts_ if false
  Vector<BaseFloat> weights_;        ///< weights (not log).
  Matrix<BaseFloat> inv_vars_;       ///< Inverted (diagonal) variances
  Matrix<BaseFloat> means_invvars_;  ///< Means times inverted variance

BaseFloat DiagGmm::ComponentLogLikelihood(const VectorBase<BaseFloat> &data,
                                          int32 comp_id) const {
  BaseFloat loglike;
  Vector<BaseFloat> data_sq(data);
  data_sq.ApplyPow(2.0);

  // loglike =  means * inv(vars) * data.
  loglike = VecVec(means_invvars_.Row(comp_id), data);
  // loglike += -0.5 * inv(vars) * data_sq.
  loglike -= 0.5 * VecVec(inv_vars_.Row(comp_id), data_sq);
  return loglike + gconsts_(comp_id);
}