plda.cc
Go to the documentation of this file.
1 // ivector/plda.cc
2 
3 // Copyright 2013 Daniel Povey
4 // 2015 David Snyder
5 
6 // See ../../COPYING for clarification regarding multiple authors
7 //
8 // Licensed under the Apache License, Version 2.0 (the "License");
9 // you may not use this file except in compliance with the License.
10 // You may obtain a copy of the License at
11 //
12 // http://www.apache.org/licenses/LICENSE-2.0
13 //
14 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
15 // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
16 // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
17 // MERCHANTABLITY OR NON-INFRINGEMENT.
18 // See the Apache 2 License for the specific language governing permissions and
19 // limitations under the License.
20 
21 #include <vector>
22 #include "ivector/plda.h"
23 
24 namespace kaldi {
25 
26 void Plda::Write(std::ostream &os, bool binary) const {
27  WriteToken(os, binary, "<Plda>");
28  mean_.Write(os, binary);
29  transform_.Write(os, binary);
30  psi_.Write(os, binary);
31  WriteToken(os, binary, "</Plda>");
32 }
33 
34 void Plda::Read(std::istream &is, bool binary) {
35  ExpectToken(is, binary, "<Plda>");
36  mean_.Read(is, binary);
37  transform_.Read(is, binary);
38  psi_.Read(is, binary);
39  ExpectToken(is, binary, "</Plda>");
41 }
42 
43 template<class Real>
47  MatrixBase<Real> *proj) {
48  int32 dim = covar.NumRows();
49  TpMatrix<Real> C(dim); // Cholesky of covar, covar = C C^T
50  C.Cholesky(covar);
51  C.Invert(); // The matrix that makes covar unit is C^{-1}, because
52  // C^{-1} covar C^{-T} = C^{-1} C C^T C^{-T} = I.
53  proj->CopyFromTp(C, kNoTrans); // set "proj" to C^{-1}.
54 }
55 
56 
58  KALDI_ASSERT(Dim() > 0);
59  offset_.Resize(Dim());
61 }
62 
63 
100  const VectorBase<double> &transformed_ivector,
101  int32 num_examples) const {
102  KALDI_ASSERT(num_examples > 0);
103  // Work out the normalization factor. The covariance for an average over
104  // "num_examples" training iVectors equals \Psi + I/num_examples.
105  Vector<double> transformed_ivector_sq(transformed_ivector);
106  transformed_ivector_sq.ApplyPow(2.0);
107  // inv_covar will equal 1.0 / (\Psi + I/num_examples).
108  Vector<double> inv_covar(psi_);
109  inv_covar.Add(1.0 / num_examples);
110  inv_covar.InvertElements();
111  // "transformed_ivector" should have covariance (\Psi + I/num_examples), i.e.
112  // within-class/num_examples plus between-class covariance. So
113  // transformed_ivector_sq . (I/num_examples + \Psi)^{-1} should be equal to
114  // the dimension.
115  double dot_prod = VecVec(inv_covar, transformed_ivector_sq);
116  return sqrt(Dim() / dot_prod);
117 }
118 
119 
120 double Plda::TransformIvector(const PldaConfig &config,
121  const VectorBase<double> &ivector,
122  int32 num_examples,
123  VectorBase<double> *transformed_ivector) const {
124  KALDI_ASSERT(ivector.Dim() == Dim() && transformed_ivector->Dim() == Dim());
125  double normalization_factor;
126  transformed_ivector->CopyFromVec(offset_);
127  transformed_ivector->AddMatVec(1.0, transform_, kNoTrans, ivector, 1.0);
128  if (config.simple_length_norm)
129  normalization_factor = sqrt(transformed_ivector->Dim())
130  / transformed_ivector->Norm(2.0);
131  else
132  normalization_factor = GetNormalizationFactor(*transformed_ivector,
133  num_examples);
134  if (config.normalize_length)
135  transformed_ivector->Scale(normalization_factor);
136  return normalization_factor;
137 }
138 
139 // "float" version of TransformIvector.
140 float Plda::TransformIvector(const PldaConfig &config,
141  const VectorBase<float> &ivector,
142  int32 num_examples,
143  VectorBase<float> *transformed_ivector) const {
144  Vector<double> tmp(ivector), tmp_out(ivector.Dim());
145  float ans = TransformIvector(config, tmp, num_examples, &tmp_out);
146  transformed_ivector->CopyFromVec(tmp_out);
147  return ans;
148 }
149 
150 
151 // There is an extended comment within this file, referencing a paper by
152 // Ioffe, that may clarify what this function is doing.
154  const VectorBase<double> &transformed_train_ivector,
155  int32 n, // number of training utterances.
156  const VectorBase<double> &transformed_test_ivector) const {
157  int32 dim = Dim();
158  double loglike_given_class, loglike_without_class;
159  { // work out loglike_given_class.
160  // "mean" will be the mean of the distribution if it comes from the
161  // training example. The mean is \frac{n \Psi}{n \Psi + I} \bar{u}^g
162  // "variance" will be the variance of that distribution, equal to
163  // I + \frac{\Psi}{n\Psi + I}.
164  Vector<double> mean(dim, kUndefined);
165  Vector<double> variance(dim, kUndefined);
166  for (int32 i = 0; i < dim; i++) {
167  mean(i) = n * psi_(i) / (n * psi_(i) + 1.0)
168  * transformed_train_ivector(i);
169  variance(i) = 1.0 + psi_(i) / (n * psi_(i) + 1.0);
170  }
171  double logdet = variance.SumLog();
172  Vector<double> sqdiff(transformed_test_ivector);
173  sqdiff.AddVec(-1.0, mean);
174  sqdiff.ApplyPow(2.0);
175  variance.InvertElements();
176  loglike_given_class = -0.5 * (logdet + M_LOG_2PI * dim +
177  VecVec(sqdiff, variance));
178  }
179  { // work out loglike_without_class. Here the mean is zero and the variance
180  // is I + \Psi.
181  Vector<double> sqdiff(transformed_test_ivector); // there is no offset.
182  sqdiff.ApplyPow(2.0);
183  Vector<double> variance(psi_);
184  variance.Add(1.0); // I + \Psi.
185  double logdet = variance.SumLog();
186  variance.InvertElements();
187  loglike_without_class = -0.5 * (logdet + M_LOG_2PI * dim +
188  VecVec(sqdiff, variance));
189  }
190  double loglike_ratio = loglike_given_class - loglike_without_class;
191  return loglike_ratio;
192 }
193 
194 
195 void Plda::SmoothWithinClassCovariance(double smoothing_factor) {
196  KALDI_ASSERT(smoothing_factor >= 0.0 && smoothing_factor <= 1.0);
197  // smoothing_factor > 1.0 is possible but wouldn't really make sense.
198 
199  KALDI_LOG << "Smoothing within-class covariance by " << smoothing_factor
200  << ", Psi is initially: " << psi_;
201  Vector<double> within_class_covar(Dim());
202  within_class_covar.Set(1.0); // It's now the current within-class covariance
203  // (a diagonal matrix) in the space transformed
204  // by transform_.
205  within_class_covar.AddVec(smoothing_factor, psi_);
210 
211  psi_.DivElements(within_class_covar);
212  KALDI_LOG << "New value of Psi is " << psi_;
213 
214  within_class_covar.ApplyPow(-0.5);
215  transform_.MulRowsVec(within_class_covar);
216 
218 }
219 
220 void Plda::ApplyTransform(const Matrix<double> &in_transform) {
221  KALDI_ASSERT(in_transform.NumRows() <= Dim()
222  && in_transform.NumCols() == Dim());
223 
224  // Apply in_transform to mean_.
225  Vector<double> mean_new(in_transform.NumRows());
226  mean_new.AddMatVec(1.0, in_transform, kNoTrans, mean_, 0.0);
227  mean_.Resize(in_transform.NumRows());
228  mean_.CopyFromVec(mean_new);
229 
230  SpMatrix<double> between_var(in_transform.NumCols()),
231  within_var(in_transform.NumCols()),
232  psi_mat(in_transform.NumCols()),
233  between_var_new(Dim()),
234  within_var_new(Dim());
235  Matrix<double> transform_invert(transform_);
236 
237  // Next, compute the between_var and within_var that existed
238  // prior to diagonalization.
239  psi_mat.AddDiagVec(1.0, psi_);
240  transform_invert.Invert();
241  within_var.AddMat2(1.0, transform_invert, kNoTrans, 0.0);
242  between_var.AddMat2Sp(1.0, transform_invert, kNoTrans, psi_mat, 0.0);
243 
244  // Next, transform the variances using the input transformation.
245  between_var_new.AddMat2Sp(1.0, in_transform, kNoTrans, between_var, 0.0);
246  within_var_new.AddMat2Sp(1.0, in_transform, kNoTrans, within_var, 0.0);
247 
248  // Finally, we need to recompute psi_ and transform_. The remainder of
249  // the code in this function is a lightly modified copy of
250  // PldaEstimator::GetOutput().
251  Matrix<double> transform1(Dim(), Dim());
252  ComputeNormalizingTransform(within_var_new, &transform1);
253  // Now transform is a matrix that if we project with it,
254  // within_var becomes unit.
255  // between_var_proj is between_var after projecting with transform1.
256  SpMatrix<double> between_var_proj(Dim());
257  between_var_proj.AddMat2Sp(1.0, transform1, kNoTrans, between_var_new, 0.0);
258 
259  Matrix<double> U(Dim(), Dim());
260  Vector<double> s(Dim());
261  // Do symmetric eigenvalue decomposition between_var_proj = U diag(s) U^T,
262  // where U is orthogonal.
263  between_var_proj.Eig(&s, &U);
264 
265  KALDI_ASSERT(s.Min() >= 0.0);
266  int32 n;
267  s.ApplyFloor(0.0, &n);
268  if (n > 0) {
269  KALDI_WARN << "Floored " << n << " eigenvalues of between-class "
270  << "variance to zero.";
271  }
272  // Sort from greatest to smallest eigenvalue.
273  SortSvd(&s, &U);
274 
275  // The transform U^T will make between_var_proj diagonal with value s
276  // (i.e. U^T U diag(s) U U^T = diag(s)). The final transform that
277  // makes within_var unit and between_var diagonal is U^T transform1,
278  // i.e. first transform1 and then U^T.
279  transform_.Resize(Dim(), Dim());
280  transform_.AddMatMat(1.0, U, kTrans, transform1, kNoTrans, 0.0);
281  psi_.Resize(Dim());
282  psi_.CopyFromVec(s);
284 }
285 
286 void PldaStats::AddSamples(double weight,
287  const Matrix<double> &group) {
288  if (dim_ == 0) {
289  Init(group.NumCols());
290  } else {
291  KALDI_ASSERT(dim_ == group.NumCols());
292  }
293  int32 n = group.NumRows(); // number of examples for this class
294  Vector<double> *mean = new Vector<double>(dim_);
295  mean->AddRowSumMat(1.0 / n, group);
296 
297  offset_scatter_.AddMat2(weight, group, kTrans, 1.0);
298  // the following statement has the same effect as if we
299  // had first subtracted the mean from each element of
300  // the group before the statement above.
301  offset_scatter_.AddVec2(-n * weight, *mean);
302 
303  class_info_.push_back(ClassInfo(weight, mean, n));
304 
305  num_classes_ ++;
306  num_examples_ += n;
307  class_weight_ += weight;
308  example_weight_ += weight * n;
309 
310  sum_.AddVec(weight, *mean);
311 }
312 
314  for (size_t i = 0; i < class_info_.size(); i++)
315  delete class_info_[i].mean;
316 }
317 
318 bool PldaStats::IsSorted() const {
319  for (size_t i = 0; i + 1 < class_info_.size(); i++)
320  if (class_info_[i+1] < class_info_[i])
321  return false;
322  return true;
323 }
324 
326  KALDI_ASSERT(dim_ == 0);
327  dim_ = dim;
328  num_classes_ = 0;
329  num_examples_ = 0;
330  class_weight_ = 0.0;
331  example_weight_ = 0.0;
332  sum_.Resize(dim);
333  offset_scatter_.Resize(dim);
334  KALDI_ASSERT(class_info_.empty());
335 }
336 
337 
339  stats_(stats) {
340  KALDI_ASSERT(stats.IsSorted());
341  InitParameters();
342 }
343 
344 
346  // Returns the part of the objf relating to offsets from the class means.
347  // within_class_count equals the sum over the classes, of the weight of that
348  // class (normally 1) times (1 - #examples) of that class, which equals the
349  // rank of the covariance we're modeling. We imagine that we're modeling (1 -
350  // #examples) separate samples, each with the within-class covariance.. the
351  // argument is a little complicated and involves an orthogonal complement of a
352  // matrix whose first row computes the mean.
353 
354  double within_class_count = stats_.example_weight_ - stats_.class_weight_,
355  within_logdet, det_sign;
356  SpMatrix<double> inv_within_var(within_var_);
357  inv_within_var.Invert(&within_logdet, &det_sign);
358  KALDI_ASSERT(det_sign == 1 && "Within-class covariance is singular");
359 
360  double objf = -0.5 * (within_class_count * (within_logdet + M_LOG_2PI * Dim())
361  + TraceSpSp(inv_within_var, stats_.offset_scatter_));
362  return objf;
363 }
364 
366  double tot_objf = 0.0;
367 
368  int32 n = -1; // the number of examples for the current class
369  SpMatrix<double> combined_inv_var(Dim());
370  // combined_inv_var = (between_var_ + within_var_ / n)^{-1}
371  double combined_var_logdet;
372 
373  for (size_t i = 0; i < stats_.class_info_.size(); i++) {
374  const ClassInfo &info = stats_.class_info_[i];
375  if (info.num_examples != n) {
376  n = info.num_examples;
377  // variance of mean of n examples is between-class + 1/n * within-class
378  combined_inv_var.CopyFromSp(between_var_);
379  combined_inv_var.AddSp(1.0 / n, within_var_);
380  combined_inv_var.Invert(&combined_var_logdet);
381  }
382  Vector<double> mean (*(info.mean));
383  mean.AddVec(-1.0 / stats_.class_weight_, stats_.sum_);
384  tot_objf += info.weight * -0.5 * (combined_var_logdet + M_LOG_2PI * Dim()
385  + VecSpVec(mean, combined_inv_var, mean));
386  }
387  return tot_objf;
388 }
389 
391  double ans1 = ComputeObjfPart1(),
392  ans2 = ComputeObjfPart2(),
393  ans = ans1 + ans2,
394  example_weights = stats_.example_weight_,
395  normalized_ans = ans / example_weights;
396  KALDI_LOG << "Within-class objf per sample is " << (ans1 / example_weights)
397  << ", between-class is " << (ans2 / example_weights)
398  << ", total is " << normalized_ans;
399  return normalized_ans;
400 }
401 
407 }
408 
411  within_var_count_ = 0.0;
413  between_var_count_ = 0.0;
414 }
415 
418  // Note: in the normal case, the expression below will be equal to the sum
419  // over the classes, of (1-n), where n is the #examples for that class. That
420  // is the rank of the scatter matrix that "offset_scatter_" has for that
421  // class. [if weights other than 1.0 are used, it will be different.]
423 }
424 
425 
471  SpMatrix<double> between_var_inv(between_var_);
472  between_var_inv.Invert();
473  SpMatrix<double> within_var_inv(within_var_);
474  within_var_inv.Invert();
475  // mixed_var will equal (between_var^{-1} + n within_var^{-1})^{-1}.
476  SpMatrix<double> mixed_var(Dim());
477  int32 n = -1; // the current number of examples for the class.
478 
479  for (size_t i = 0; i < stats_.class_info_.size(); i++) {
480  const ClassInfo &info = stats_.class_info_[i];
481  double weight = info.weight;
482  if (info.num_examples != n) {
483  n = info.num_examples;
484  mixed_var.CopyFromSp(between_var_inv);
485  mixed_var.AddSp(n, within_var_inv);
486  mixed_var.Invert();
487  }
488  Vector<double> m = *(info.mean); // the mean for this class.
489  m.AddVec(-1.0 / stats_.class_weight_, stats_.sum_); // remove global mean
490  Vector<double> temp(Dim()); // n within_var^{-1} m
491  temp.AddSpVec(n, within_var_inv, m, 0.0);
492  Vector<double> w(Dim()); // w, as defined in the comment.
493  w.AddSpVec(1.0, mixed_var, temp, 0.0);
494  Vector<double> m_w(m); // m - w
495  m_w.AddVec(-1.0, w);
496  between_var_stats_.AddSp(weight, mixed_var);
497  between_var_stats_.AddVec2(weight, w);
498  between_var_count_ += weight;
499  within_var_stats_.AddSp(weight * n, mixed_var);
500  within_var_stats_.AddVec2(weight * n, m_w);
501  within_var_count_ += weight;
502  }
503 }
504 
510 
511  KALDI_LOG << "Trace of within-class variance is " << within_var_.Trace();
512  KALDI_LOG << "Trace of between-class variance is " << between_var_.Trace();
513 }
514 
515 
521  KALDI_VLOG(2) << "Objective function is " << ComputeObjf();
522 }
523 
524 
526  Plda *plda) {
527  KALDI_ASSERT(stats_.example_weight_ > 0 && "Cannot estimate with no stats");
528  for (int32 i = 0; i < config.num_em_iters; i++) {
529  KALDI_LOG << "Plda estimation iteration " << i
530  << " of " << config.num_em_iters;
531  EstimateOneIter();
532  }
533  GetOutput(plda);
534 }
535 
536 
538  plda->mean_ = stats_.sum_;
539  plda->mean_.Scale(1.0 / stats_.class_weight_);
540  KALDI_LOG << "Norm of mean of iVector distribution is "
541  << plda->mean_.Norm(2.0);
542 
543  Matrix<double> transform1(Dim(), Dim());
545  // now transform is a matrix that if we project with it,
546  // within_var_ becomes unit.
547 
548  // between_var_proj is between_var after projecting with transform1.
549  SpMatrix<double> between_var_proj(Dim());
550  between_var_proj.AddMat2Sp(1.0, transform1, kNoTrans, between_var_, 0.0);
551 
552  Matrix<double> U(Dim(), Dim());
553  Vector<double> s(Dim());
554  // Do symmetric eigenvalue decomposition between_var_proj = U diag(s) U^T,
555  // where U is orthogonal.
556  between_var_proj.Eig(&s, &U);
557 
558  KALDI_ASSERT(s.Min() >= 0.0);
559  int32 n;
560  s.ApplyFloor(0.0, &n);
561  if (n > 0) {
562  KALDI_WARN << "Floored " << n << " eigenvalues of between-class "
563  << "variance to zero.";
564  }
565  // Sort from greatest to smallest eigenvalue.
566  SortSvd(&s, &U);
567 
568  // The transform U^T will make between_var_proj diagonal with value s
569  // (i.e. U^T U diag(s) U U^T = diag(s)). The final transform that
570  // makes within_var_ unit and between_var_ diagonal is U^T transform1,
571  // i.e. first transform1 and then U^T.
572 
573  plda->transform_.Resize(Dim(), Dim());
574  plda->transform_.AddMatMat(1.0, U, kTrans, transform1, kNoTrans, 0.0);
575  plda->psi_ = s;
576 
577  KALDI_LOG << "Diagonal of between-class variance in normalized space is " << s;
578 
579  if (GetVerboseLevel() >= 2) { // at higher verbose levels, do a self-test
580  // (just tests that this function does what it
581  // should).
582  SpMatrix<double> tmp_within(Dim());
583  tmp_within.AddMat2Sp(1.0, plda->transform_, kNoTrans, within_var_, 0.0);
584  KALDI_ASSERT(tmp_within.IsUnit(0.0001));
585  SpMatrix<double> tmp_between(Dim());
586  tmp_between.AddMat2Sp(1.0, plda->transform_, kNoTrans, between_var_, 0.0);
587  KALDI_ASSERT(tmp_between.IsDiagonal(0.0001));
588  Vector<double> psi(Dim());
589  psi.CopyDiagFromSp(tmp_between);
590  AssertEqual(psi, plda->psi_);
591  }
592  plda->ComputeDerivedVars();
593 }
594 
596  const Vector<double> &ivector) {
597  if (mean_stats_.Dim() == 0) {
598  mean_stats_.Resize(ivector.Dim());
599  variance_stats_.Resize(ivector.Dim());
600  }
601  KALDI_ASSERT(weight >= 0.0);
602  tot_weight_ += weight;
603  mean_stats_.AddVec(weight, ivector);
604  variance_stats_.AddVec2(weight, ivector);
605 }
606 
608  const Vector<float> &ivector) {
609  Vector<double> ivector_dbl(ivector);
610  this->AddStats(weight, ivector_dbl);
611 }
612 
614  Plda *plda) const {
615  KALDI_ASSERT(tot_weight_ > 0.0);
616  int32 dim = mean_stats_.Dim();
617  KALDI_ASSERT(dim == plda->Dim());
618  Vector<double> mean(mean_stats_);
619  mean.Scale(1.0 / tot_weight_);
620  SpMatrix<double> variance(variance_stats_);
621  variance.Scale(1.0 / tot_weight_);
622  variance.AddVec2(-1.0, mean); // Make it the uncentered variance.
623 
624  // mean_diff of the adaptation data from the training data. We optionally add
625  // this to our total covariance matrix
626  Vector<double> mean_diff(mean);
627  mean_diff.AddVec(-1.0, plda->mean_);
628  KALDI_ASSERT(config.mean_diff_scale >= 0.0);
629  variance.AddVec2(config.mean_diff_scale, mean_diff);
630 
631  // update the plda's mean data-member with our adaptation-data mean.
632  plda->mean_.CopyFromVec(mean);
633 
634 
635  // transform_model_ is a row-scaled version of plda->transform_ that
636  // transforms into the space where the total covariance is 1.0. Because
637  // plda->transform_ transforms into a space where the within-class covar is
638  // 1.0 and the the between-class covar is diag(plda->psi_), we need to scale
639  // each dimension i by 1.0 / sqrt(1.0 + plda->psi_(i))
640 
641  Matrix<double> transform_mod(plda->transform_);
642  for (int32 i = 0; i < dim; i++)
643  transform_mod.Row(i).Scale(1.0 / sqrt(1.0 + plda->psi_(i)));
644 
645  // project the variance of the adaptation set into this space where
646  // the total covariance is unit.
647  SpMatrix<double> variance_proj(dim);
648  variance_proj.AddMat2Sp(1.0, transform_mod, kNoTrans,
649  variance, 0.0);
650 
651  // Do eigenvalue decomposition of variance_proj; this will tell us the
652  // directions in which the adaptation-data covariance is more than
653  // the training-data covariance.
654  Matrix<double> P(dim, dim);
655  Vector<double> s(dim);
656  variance_proj.Eig(&s, &P);
657  SortSvd(&s, &P);
658  KALDI_LOG << "Eigenvalues of adaptation-data total-covariance in space where "
659  << "training-data total-covariance is unit, is: " << s;
660 
661  // W, B are the (within,between)-class covars in the space transformed by
662  // transform_mod.
663  SpMatrix<double> W(dim), B(dim);
664  for (int32 i = 0; i < dim; i++) {
665  W(i, i) = 1.0 / (1.0 + plda->psi_(i)),
666  B(i, i) = plda->psi_(i) / (1.0 + plda->psi_(i));
667  }
668 
669  // OK, so variance_proj (projected by transform_mod) is P diag(s) P^T.
670  // Suppose that after transform_mod we project by P^T. Then the adaptation-data's
671  // variance would be P^T P diag(s) P^T P = diag(s), and the PLDA model's
672  // within class variance would be P^T W P and its between-class variance would be
673  // P^T B P. We'd still have that W+B = I in this space.
674  // First let's compute these projected variances... we call the "proj2" because
675  // it's after the data has been projected twice (actually, transformed, as there is no
676  // dimension loss), by transform_mod and then P^T.
677 
678  SpMatrix<double> Wproj2(dim), Bproj2(dim);
679  Wproj2.AddMat2Sp(1.0, P, kTrans, W, 0.0);
680  Bproj2.AddMat2Sp(1.0, P, kTrans, B, 0.0);
681 
682  Matrix<double> Ptrans(P, kTrans);
683 
684  SpMatrix<double> Wproj2mod(Wproj2), Bproj2mod(Bproj2);
685 
686  for (int32 i = 0; i < dim; i++) {
687  // For this eigenvalue, compute the within-class covar projected with this direction,
688  // and the same for between.
689  BaseFloat within = Wproj2(i, i),
690  between = Bproj2(i, i);
691  KALDI_LOG << "For " << i << "'th eigenvalue, value is " << s(i)
692  << ", within-class covar in this direction is " << within
693  << ", between-class is " << between;
694  if (s(i) > 1.0) {
695  double excess_eig = s(i) - 1.0;
696  double excess_within_covar = excess_eig * config.within_covar_scale,
697  excess_between_covar = excess_eig * config.between_covar_scale;
698  Wproj2mod(i, i) += excess_within_covar;
699  Bproj2mod(i, i) += excess_between_covar;
700  } /*
701  Below I was considering a method like below, to try to scale up
702  the dimensions that had less variance than expected in our sample..
703  this didn't help, and actually when I set that power to +0.2 instead
704  of -0.5 it gave me an improvement on sre08. But I'm not sure
705  about this.. it just doesn't seem right.
706  else {
707  BaseFloat scale = pow(std::max(1.0e-10, s(i)), -0.5);
708  BaseFloat max_scale = 10.0; // I'll make this configurable later.
709  scale = std::min(scale, max_scale);
710  Ptrans.Row(i).Scale(scale);
711  } */
712  }
713 
714  // combined transform "transform_mod" and then P^T that takes us to the space
715  // where {W,B}proj2{,mod} are.
716  Matrix<double> combined_trans(dim, dim);
717  combined_trans.AddMatMat(1.0, Ptrans, kNoTrans,
718  transform_mod, kNoTrans, 0.0);
719  Matrix<double> combined_trans_inv(combined_trans); // ... and its inverse.
720  combined_trans_inv.Invert();
721 
722  // Wmod and Bmod are as Wproj2 and Bproj2 but taken back into the original
723  // iVector space.
724  SpMatrix<double> Wmod(dim), Bmod(dim);
725  Wmod.AddMat2Sp(1.0, combined_trans_inv, kNoTrans, Wproj2mod, 0.0);
726  Bmod.AddMat2Sp(1.0, combined_trans_inv, kNoTrans, Bproj2mod, 0.0);
727 
728  TpMatrix<double> C(dim);
729  // Do Cholesky Wmod = C C^T. Now if we use C^{-1} as a transform, we have
730  // C^{-1} W C^{-T} = I, so it makes the within-class covar unit.
731  C.Cholesky(Wmod);
732  TpMatrix<double> Cinv(C);
733  Cinv.Invert();
734 
735  // Bmod_proj is Bmod projected by Cinv.
736  SpMatrix<double> Bmod_proj(dim);
737  Bmod_proj.AddTp2Sp(1.0, Cinv, kNoTrans, Bmod, 0.0);
738  Vector<double> psi_new(dim);
739  Matrix<double> Q(dim, dim);
740  // Do symmetric eigenvalue decomposition of Bmod_proj, so
741  // Bmod_proj = Q diag(psi_new) Q^T
742  Bmod_proj.Eig(&psi_new, &Q);
743  SortSvd(&psi_new, &Q);
744  // This means that if we use Q^T as a transform, then Q^T Bmod_proj Q =
745  // diag(psi_new), hence Q^T diagonalizes Bmod_proj (while leaving the
746  // within-covar unit).
747  // The final transform we want, that projects from our original
748  // space to our newly normalized space, is:
749  // first Cinv, then Q^T, i.e. the
750  // matrix Q^T Cinv.
751  Matrix<double> final_transform(dim, dim);
752  final_transform.AddMatTp(1.0, Q, kTrans, Cinv, kNoTrans, 0.0);
753 
754  KALDI_LOG << "Old diagonal of between-class covar was: "
755  << plda->psi_ << ", new diagonal is "
756  << psi_new;
757  plda->transform_.CopyFromMat(final_transform);
758  plda->psi_.CopyFromVec(psi_new);
759 }
760 
761 } // namespace kaldi
This code computes Goodness of Pronunciation (GOP) and extracts phone-level pronunciation feature for...
Definition: chain.dox:20
bool IsUnit(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:480
void Estimate(const PldaEstimationConfig &config, Plda *output)
Definition: plda.cc:525
Packed symetric matrix class.
Definition: matrix-common.h:62
SpMatrix< double > between_var_
Definition: plda.h:278
void Write(std::ostream &out, bool binary) const
write to stream.
#define M_LOG_2PI
Definition: kaldi-math.h:60
void Scale(Real c)
void AddRowSumMat(Real alpha, const MatrixBase< Real > &M, Real beta=1.0)
Does *this = alpha * (sum of rows of M) + beta * *this.
SpMatrix< double > within_var_stats_
Definition: plda.h:281
MatrixIndexT NumCols() const
Returns number of columns (or zero for empty matrix).
Definition: kaldi-matrix.h:67
Base class which provides matrix operations not involving resizing or allocation. ...
Definition: kaldi-matrix.h:49
double example_weight_
Definition: plda.h:199
Real Trace() const
Definition: sp-matrix.cc:171
double TransformIvector(const PldaConfig &config, const VectorBase< double > &ivector, int32 num_enroll_examples, VectorBase< double > *transformed_ivector) const
Transforms an iVector into a space where the within-class variance is unit and between-class variance...
Definition: plda.cc:120
int32 GetVerboseLevel()
Get verbosity level, usually set via command line &#39;–verbose=&#39; switch.
Definition: kaldi-error.h:60
void Read(std::istream &is, bool binary)
Definition: plda.cc:34
void Write(std::ostream &Out, bool binary) const
Writes to C++ stream (option to write in binary).
Real SumLog() const
Returns sum of the logs of the elements.
Vector< double > offset_
Definition: plda.h:155
Matrix< double > transform_
Definition: plda.h:149
kaldi::int32 int32
void Eig(VectorBase< Real > *s, MatrixBase< Real > *P=NULL) const
Solves the symmetric eigenvalue problem: at end we should have (*this) = P * diag(s) * P^T...
Definition: qr.cc:433
SpMatrix< double > offset_scatter_
Definition: plda.h:204
Vector< double > sum_
Definition: plda.h:201
void SetUnit()
< Set to zero
void Resize(MatrixIndexT length, MatrixResizeType resize_type=kSetZero)
Set vector to a specified size (can be zero).
void Init(int32 dim)
Definition: plda.cc:325
void AddSpVec(const Real alpha, const SpMatrix< Real > &M, const VectorBase< Real > &v, const Real beta)
Add symmetric positive definite matrix times vector: this <– beta*this + alpha*M*v.
const PldaStats & stats_
Definition: plda.h:275
void CopyFromMat(const MatrixBase< OtherReal > &M, MatrixTransposeType trans=kNoTrans)
Copy given matrix. (no resize is done).
Real Min() const
Returns the minimum value of any element, or +infinity for the empty vector.
MatrixIndexT NumRows() const
void CopyFromSp(const SpMatrix< Real > &other)
Definition: sp-matrix.h:85
Real Norm(Real p) const
Compute the p-th norm of the vector.
void AddTp2Sp(const Real alpha, const TpMatrix< Real > &T, MatrixTransposeType transM, const SpMatrix< Real > &A, const Real beta=0.0)
The following function does: this <– beta*this + alpha * T * A * T^T.
Definition: sp-matrix.cc:1139
static void ComputeNormalizingTransform(const SpMatrix< Real > &covar, MatrixBase< Real > *proj)
This function computes a projection matrix that when applied makes the covariance unit (i...
Definition: plda.cc:46
void EstimateFromStats()
Definition: plda.cc:505
void ApplyFloor(Real floor_val, MatrixIndexT *floored_count=nullptr)
Applies floor to all elements.
Definition: kaldi-vector.h:149
void CopyFromVec(const VectorBase< Real > &v)
Copy data from another vector (must match own size).
void AddVec2(const Real alpha, const VectorBase< OtherReal > &v)
rank-one update, this <– this + alpha v v&#39;
Definition: sp-matrix.cc:946
void Read(std::istream &in, bool binary, bool add=false)
read from stream.
void Cholesky(const SpMatrix< Real > &orig)
Definition: tp-matrix.cc:88
void ResetPerIterStats()
Definition: plda.cc:409
const SubVector< Real > Row(MatrixIndexT i) const
Return specific row of matrix [const].
Definition: kaldi-matrix.h:188
Vector< double > psi_
Definition: plda.h:152
PldaEstimator(const PldaStats &stats)
Definition: plda.cc:338
void AddSp(const Real alpha, const SpMatrix< Real > &Ma)
Definition: sp-matrix.h:211
void ExpectToken(std::istream &is, bool binary, const char *token)
ExpectToken tries to read in the given token, and throws an exception on failure. ...
Definition: io-funcs.cc:191
double class_weight_
Definition: plda.h:198
struct rnnlm::@11::@12 n
void SmoothWithinClassCovariance(double smoothing_factor)
This function smooths the within-class covariance by adding to it, smoothing_factor (e...
Definition: plda.cc:195
void CopyFromTp(const TpMatrix< OtherReal > &M, MatrixTransposeType trans=kNoTrans)
Copy given tpmatrix. (no resize is done).
void EstimateOneIter()
Definition: plda.cc:516
void AddMatMat(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const Real beta)
double TraceSpSp(const SpMatrix< double > &A, const SpMatrix< double > &B)
Definition: sp-matrix.cc:326
Real VecSpVec(const VectorBase< Real > &v1, const SpMatrix< Real > &M, const VectorBase< Real > &v2)
Computes v1^T * M * v2.
Definition: sp-matrix.cc:964
Packed symetric matrix class.
Definition: matrix-common.h:63
#define KALDI_WARN
Definition: kaldi-error.h:150
void Write(std::ostream &os, bool binary) const
Definition: plda.cc:26
double LogLikelihoodRatio(const VectorBase< double > &transformed_enroll_ivector, int32 num_enroll_utts, const VectorBase< double > &transformed_test_ivector) const
Returns the log-likelihood ratio log (p(test_ivector | same) / p(test_ivector | different)).
Definition: plda.cc:153
void WriteToken(std::ostream &os, bool binary, const char *token)
The WriteToken functions are for writing nonempty sequences of non-space characters.
Definition: io-funcs.cc:134
MatrixIndexT Dim() const
Returns the dimension of the vector.
Definition: kaldi-vector.h:64
Vector< double > * mean
Definition: plda.h:210
void GetStatsFromClassMeans()
GetStatsFromClassMeans() is the more complicated part of PLDA estimation.
Definition: plda.cc:470
void Scale(Real alpha)
Multiplies all elements by this constant.
void AddMatVec(const Real alpha, const MatrixBase< Real > &M, const MatrixTransposeType trans, const VectorBase< Real > &v, const Real beta)
Add matrix times vector : this <– beta*this + alpha*M*v.
Definition: kaldi-vector.cc:92
void ComputeDerivedVars()
Definition: plda.cc:57
double GetNormalizationFactor(const VectorBase< double > &transformed_ivector, int32 num_examples) const
This returns a normalization factor, which is a quantity we must multiply "transformed_ivector" by so...
Definition: plda.cc:99
double ComputeObjfPart2() const
Returns the part of the obj relating to the class means (total_not normalized)
Definition: plda.cc:365
void AddStats(double weight, const Vector< double > &ivector)
Definition: plda.cc:595
void MulRowsVec(const VectorBase< Real > &scale)
Equivalent to (*this) = diag(scale) * (*this).
double ComputeObjfPart1() const
Returns the part of the objf relating to offsets from the class means.
Definition: plda.cc:345
void GetOutput(Plda *plda)
Definition: plda.cc:537
void UpdatePlda(const PldaUnsupervisedAdaptorConfig &config, Plda *plda) const
Definition: plda.cc:613
void InvertElements()
Invert all elements.
void InitParameters()
Definition: plda.cc:402
double between_var_count_
Definition: plda.h:284
Vector< double > mean_
Definition: plda.h:148
void ApplyTransform(const Matrix< double > &in_transform)
Apply a transform to the PLDA model.
Definition: plda.cc:220
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
int32 Dim() const
Definition: plda.h:255
void AddSamples(double weight, const Matrix< double > &group)
The dimension is set up the first time you add samples.
Definition: plda.cc:286
double within_var_count_
Definition: plda.h:282
SpMatrix< double > within_var_
Definition: plda.h:277
MatrixIndexT NumRows() const
Returns number of rows (or zero for empty matrix).
Definition: kaldi-matrix.h:64
void Set(Real f)
Set all members of a vector to a specified value.
void AddMat2Sp(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const SpMatrix< Real > &A, const Real beta=0.0)
Extension of rank-N update: this <– beta*this + alpha * M * A * M^T.
Definition: sp-matrix.cc:982
void ApplyPow(Real power)
Take all elements of vector to a power.
Definition: kaldi-vector.h:179
std::vector< ClassInfo > class_info_
Definition: plda.h:220
static void AssertEqual(float a, float b, float relative_tolerance=0.001)
assert abs(a - b) <= relative_tolerance * (abs(a)+abs(b))
Definition: kaldi-math.h:276
#define KALDI_VLOG(v)
Definition: kaldi-error.h:156
void Resize(const MatrixIndexT r, const MatrixIndexT c, MatrixResizeType resize_type=kSetZero, MatrixStrideType stride_type=kDefaultStride)
Sets matrix to a specified size (zero is OK as long as both r and c are zero).
void Resize(MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
Definition: sp-matrix.h:81
bool IsSorted() const
Definition: plda.cc:318
SpMatrix< double > between_var_stats_
Definition: plda.h:283
void Invert(Real *log_det=NULL, Real *det_sign=NULL, bool inverse_needed=true)
matrix inverse.
Definition: kaldi-matrix.cc:38
void AddMatTp(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const TpMatrix< Real > &B, MatrixTransposeType transB, const Real beta)
this <– beta*this + alpha*A*B.
Definition: kaldi-matrix.h:725
Provides a vector abstraction class.
Definition: kaldi-vector.h:41
void Add(Real c)
Add a constant to each element of a vector.
int32 Dim() const
Definition: plda.h:140
bool normalize_length
Definition: plda.h:53
void Invert(Real *logdet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
matrix inverse.
Definition: sp-matrix.cc:219
#define KALDI_LOG
Definition: kaldi-error.h:153
Real VecVec(const VectorBase< Real > &a, const VectorBase< Real > &b)
Returns dot product between v1 and v2.
Definition: kaldi-vector.cc:37
void AddVec(const Real alpha, const VectorBase< OtherReal > &v)
Add vector : *this = *this + alpha * rv (with casting between floats and doubles) ...
void Read(std::istream &in, bool binary, bool add=false)
Read function using C++ streams.
double ComputeObjf() const
Returns the objective-function per sample.
Definition: plda.cc:390
void SortSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt, bool sort_on_absolute_value)
Function to ensure that SVD is sorted.
bool simple_length_norm
Definition: plda.h:54
void GetStatsFromIntraClass()
Definition: plda.cc:416