sp-matrix.cc
Go to the documentation of this file.
1 // matrix/sp-matrix.cc
2 
3 // Copyright 2009-2011 Lukas Burget; Ondrej Glembek; Microsoft Corporation
4 // Saarland University; Petr Schwarz; Yanmin Qian;
5 // Haihua Xu
6 
7 // See ../../COPYING for clarification regarding multiple authors
8 //
9 // Licensed under the Apache License, Version 2.0 (the "License");
10 // you may not use this file except in compliance with the License.
11 // You may obtain a copy of the License at
12 
13 // http://www.apache.org/licenses/LICENSE-2.0
14 
15 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16 // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
17 // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
18 // MERCHANTABLITY OR NON-INFRINGEMENT.
19 // See the Apache 2 License for the specific language governing permissions and
20 // limitations under the License.
21 
22 #include <limits>
23 
24 #include "matrix/sp-matrix.h"
25 #include "matrix/kaldi-vector.h"
26 #include "matrix/kaldi-matrix.h"
28 #include "matrix/cblas-wrappers.h"
29 
30 namespace kaldi {
31 
32 // ****************************************************************************
33 // Returns the log-determinant if +ve definite, else KALDI_ERR.
34 // ****************************************************************************
35 template<typename Real>
37  TpMatrix<Real> chol(this->NumRows());
38  double det = 0.0;
39  double diag;
40  chol.Cholesky(*this); // Will throw exception if not +ve definite!
41 
42  for (MatrixIndexT i = 0; i < this->NumRows(); i++) {
43  diag = static_cast<double>(chol(i, i));
44  det += kaldi::Log(diag);
45  }
46  return static_cast<Real>(2*det);
47 }
48 
49 
50 template<typename Real>
52  std::swap(this->data_, other->data_);
53  std::swap(this->num_rows_, other->num_rows_);
54 }
55 
56 template<typename Real>
59  Real tolerance) const {
60  Eig(s, P);
61  Real max = s->Max(), min = s->Min();
62  KALDI_ASSERT(-min <= tolerance * max);
63  s->ApplyFloor(0.0);
64 }
65 
66 template<typename Real>
68  Vector<Real> s(this->NumRows());
69  this->Eig(&s, static_cast<MatrixBase<Real>*>(NULL));
70  return std::max(s.Max(), -s.Min());
71 }
72 
73 // returns true if positive definite--uses cholesky.
74 template<typename Real>
76  MatrixIndexT D = (*this).NumRows();
77  KALDI_ASSERT(D > 0);
78  try {
79  TpMatrix<Real> C(D);
80  C.Cholesky(*this);
81  for (MatrixIndexT r = 0; r < D; r++)
82  if (C(r, r) == 0.0) return false;
83  return true;
84  }
85  catch(...) { // not positive semidefinite.
86  return false;
87  }
88 }
89 
90 template<typename Real>
91 void SpMatrix<Real>::ApplyPow(Real power) {
92  if (power == 1) return; // can do nothing.
93  MatrixIndexT D = this->NumRows();
94  KALDI_ASSERT(D > 0);
95  Matrix<Real> U(D, D);
96  Vector<Real> l(D);
97  (*this).SymPosSemiDefEig(&l, &U);
98 
99  Vector<Real> l_copy(l);
100  try {
101  l.ApplyPow(power * 0.5);
102  }
103  catch(...) {
104  KALDI_ERR << "Error taking power " << (power * 0.5) << " of vector "
105  << l_copy;
106  }
107  U.MulColsVec(l);
108  (*this).AddMat2(1.0, U, kNoTrans, 0.0);
109 }
110 
111 template<typename Real>
113  SpCopyType copy_type) {
114  KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
115  MatrixIndexT D = this->NumRows();
116 
117  switch (copy_type) {
118  case kTakeMeanAndCheck:
119  {
120  Real good_sum = 0.0, bad_sum = 0.0;
121  for (MatrixIndexT i = 0; i < D; i++) {
122  for (MatrixIndexT j = 0; j < i; j++) {
123  Real a = M(i, j), b = M(j, i), avg = 0.5*(a+b), diff = 0.5*(a-b);
124  (*this)(i, j) = avg;
125  good_sum += std::abs(avg);
126  bad_sum += std::abs(diff);
127  }
128  good_sum += std::abs(M(i, i));
129  (*this)(i, i) = M(i, i);
130  }
131  if (bad_sum > 0.01 * good_sum) {
132  KALDI_ERR << "SpMatrix::Copy(), source matrix is not symmetric: "
133  << bad_sum << ">" << good_sum;
134  }
135  break;
136  }
137  case kTakeMean:
138  {
139  for (MatrixIndexT i = 0; i < D; i++) {
140  for (MatrixIndexT j = 0; j < i; j++) {
141  (*this)(i, j) = 0.5*(M(i, j) + M(j, i));
142  }
143  (*this)(i, i) = M(i, i);
144  }
145  break;
146  }
147  case kTakeLower:
148  { // making this one a bit more efficient.
149  const Real *src = M.Data();
150  Real *dest = this->data_;
151  MatrixIndexT stride = M.Stride();
152  for (MatrixIndexT i = 0; i < D; i++) {
153  for (MatrixIndexT j = 0; j <= i; j++)
154  dest[j] = src[j];
155  dest += i + 1;
156  src += stride;
157  }
158  }
159  break;
160  case kTakeUpper:
161  for (MatrixIndexT i = 0; i < D; i++)
162  for (MatrixIndexT j = 0; j <= i; j++)
163  (*this)(i, j) = M(j, i);
164  break;
165  default:
166  KALDI_ASSERT("Invalid argument to SpMatrix::CopyFromMat");
167  }
168 }
169 
170 template<typename Real>
171 Real SpMatrix<Real>::Trace() const {
172  const Real *data = this->data_;
173  MatrixIndexT num_rows = this->num_rows_;
174  Real ans = 0.0;
175  for (int32 i = 1; i <= num_rows; i++, data += i)
176  ans += *data;
177  return ans;
178 }
179 
180 // diagonal update, this <-- this + diag(v)
181 template<typename Real>
182 template<typename OtherReal>
183 void SpMatrix<Real>::AddDiagVec(const Real alpha, const VectorBase<OtherReal> &v) {
184  int32 num_rows = this->num_rows_;
185  KALDI_ASSERT(num_rows == v.Dim() && num_rows > 0);
186  const OtherReal *src = v.Data();
187  Real *dst = this->data_;
188  if (alpha == 1.0)
189  for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
190  *dst += *src;
191  else
192  for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
193  *dst += alpha * *src;
194 }
195 
196 // instantiate the template above.
197 template
198 void SpMatrix<float>::AddDiagVec(const float alpha,
199  const VectorBase<double> &v);
200 
201 template
202 void SpMatrix<double>::AddDiagVec(const double alpha,
203  const VectorBase<float> &v);
204 
205 template
206 void SpMatrix<float>::AddDiagVec(const float alpha,
207  const VectorBase<float> &v);
208 
209 template
210 void SpMatrix<double>::AddDiagVec(const double alpha,
211  const VectorBase<double> &v);
212 
213 template<>
214 template<>
215 void SpMatrix<double>::AddVec2(const double alpha, const VectorBase<double> &v);
216 
217 #ifndef HAVE_ATLAS
218 template<typename Real>
219 void SpMatrix<Real>::Invert(Real *logdet, Real *det_sign, bool need_inverse) {
220  // these are CLAPACK types
221  KaldiBlasInt result;
222  KaldiBlasInt rows = static_cast<int>(this->num_rows_);
223  KaldiBlasInt* p_ipiv = new KaldiBlasInt[rows];
224  Real *p_work; // workspace for the lapack function
225  void *temp;
226  if ((p_work = static_cast<Real*>(
227  KALDI_MEMALIGN(16, sizeof(Real) * rows, &temp))) == NULL) {
228  delete[] p_ipiv;
229  throw std::bad_alloc();
230  }
231 #ifdef HAVE_OPENBLAS
232  memset(p_work, 0, sizeof(Real) * rows); // gets rid of a probably
233  // spurious Valgrind warning about jumps depending upon uninitialized values.
234 #endif
235 
236 
237  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
238  // of the data. We have a row-wise storage, therefore, we need to "invert"
239  clapack_Xsptrf(&rows, this->data_, p_ipiv, &result);
240 
241 
242  KALDI_ASSERT(result >= 0 && "Call to CLAPACK ssptrf_ called with wrong arguments");
243 
244  if (result > 0) { // Singular...
245  if (det_sign) *det_sign = 0;
246  if (logdet) *logdet = -std::numeric_limits<Real>::infinity();
247  if (need_inverse) KALDI_ERR << "CLAPACK stptrf_ : factorization failed";
248  } else { // Not singular.. compute log-determinant if needed.
249  if (logdet != NULL || det_sign != NULL) {
250  Real prod = 1.0, log_prod = 0.0;
251  int sign = 1;
252  for (int i = 0; i < (int)this->num_rows_; i++) {
253  if (p_ipiv[i] > 0) { // not a 2x2 block...
254  // if (p_ipiv[i] != i+1) sign *= -1; // row swap.
255  Real diag = (*this)(i, i);
256  prod *= diag;
257  } else { // negative: 2x2 block. [we are in first of the two].
258  i++; // skip over the first of the pair.
259  // each 2x2 block...
260  Real diag1 = (*this)(i, i), diag2 = (*this)(i-1, i-1),
261  offdiag = (*this)(i, i-1);
262  Real thisdet = diag1*diag2 - offdiag*offdiag;
263  // thisdet == determinant of 2x2 block.
264  // The following line is more complex than it looks: there are 2 offsets of
265  // 1 that cancel.
266  prod *= thisdet;
267  }
268  if (i == (int)(this->num_rows_-1) || fabs(prod) < 1.0e-10 || fabs(prod) > 1.0e+10) {
269  if (prod < 0) { prod = -prod; sign *= -1; }
270  log_prod += kaldi::Log(std::abs(prod));
271  prod = 1.0;
272  }
273  }
274  if (logdet != NULL) *logdet = log_prod;
275  if (det_sign != NULL) *det_sign = sign;
276  }
277  }
278  if (!need_inverse) {
279  delete [] p_ipiv;
280  KALDI_MEMALIGN_FREE(p_work);
281  return; // Don't need what is computed next.
282  }
283  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
284  // of the data. We have a row-wise storage, therefore, we need to "invert"
285  clapack_Xsptri(&rows, this->data_, p_ipiv, p_work, &result);
286 
287  KALDI_ASSERT(result >=0 &&
288  "Call to CLAPACK ssptri_ called with wrong arguments");
289 
290  if (result != 0) {
291  KALDI_ERR << "CLAPACK ssptrf_ : Matrix is singular";
292  }
293 
294  delete [] p_ipiv;
295  KALDI_MEMALIGN_FREE(p_work);
296 }
297 #else
298 // in the ATLAS case, these are not implemented using a library and we back off to something else.
299 template<typename Real>
300 void SpMatrix<Real>::Invert(Real *logdet, Real *det_sign, bool need_inverse) {
301  Matrix<Real> M(this->NumRows(), this->NumCols());
302  M.CopyFromSp(*this);
303  M.Invert(logdet, det_sign, need_inverse);
304  if (need_inverse)
305  for (MatrixIndexT i = 0; i < this->NumRows(); i++)
306  for (MatrixIndexT j = 0; j <= i; j++)
307  (*this)(i, j) = M(i, j);
308 }
309 #endif
310 
311 template<typename Real>
312 void SpMatrix<Real>::InvertDouble(Real *logdet, Real *det_sign,
313  bool inverse_needed) {
314  SpMatrix<double> dmat(*this);
315  double logdet_tmp, det_sign_tmp;
316  dmat.Invert(logdet ? &logdet_tmp : NULL,
317  det_sign ? &det_sign_tmp : NULL,
318  inverse_needed);
319  if (logdet) *logdet = logdet_tmp;
320  if (det_sign) *det_sign = det_sign_tmp;
321  (*this).CopyFromSp(dmat);
322 }
323 
324 
325 
326 double TraceSpSp(const SpMatrix<double> &A, const SpMatrix<double> &B) {
327  KALDI_ASSERT(A.NumRows() == B.NumRows());
328  const double *Aptr = A.Data();
329  const double *Bptr = B.Data();
330  MatrixIndexT R = A.NumRows();
331  MatrixIndexT RR = (R * (R + 1)) / 2;
332  double all_twice = 2.0 * cblas_Xdot(RR, Aptr, 1, Bptr, 1);
333  // "all_twice" contains twice the vector-wise dot-product... this is
334  // what we want except the diagonal elements are represented
335  // twice.
336  double diag_once = 0.0;
337  for (MatrixIndexT row_plus_two = 2; row_plus_two <= R + 1; row_plus_two++) {
338  diag_once += *Aptr * *Bptr;
339  Aptr += row_plus_two;
340  Bptr += row_plus_two;
341  }
342  return all_twice - diag_once;
343 }
344 
345 
346 float TraceSpSp(const SpMatrix<float> &A, const SpMatrix<float> &B) {
347  KALDI_ASSERT(A.NumRows() == B.NumRows());
348  const float *Aptr = A.Data();
349  const float *Bptr = B.Data();
350  MatrixIndexT R = A.NumRows();
351  MatrixIndexT RR = (R * (R + 1)) / 2;
352  float all_twice = 2.0 * cblas_Xdot(RR, Aptr, 1, Bptr, 1);
353  // "all_twice" contains twice the vector-wise dot-product... this is
354  // what we want except the diagonal elements are represented
355  // twice.
356  float diag_once = 0.0;
357  for (MatrixIndexT row_plus_two = 2; row_plus_two <= R + 1; row_plus_two++) {
358  diag_once += *Aptr * *Bptr;
359  Aptr += row_plus_two;
360  Bptr += row_plus_two;
361  }
362  return all_twice - diag_once;
363 }
364 
365 
366 template<typename Real, typename OtherReal>
367 Real TraceSpSp(const SpMatrix<Real> &A, const SpMatrix<OtherReal> &B) {
368  KALDI_ASSERT(A.NumRows() == B.NumRows());
369  Real ans = 0.0;
370  const Real *Aptr = A.Data();
371  const OtherReal *Bptr = B.Data();
372  MatrixIndexT row, col, R = A.NumRows();
373  for (row = 0; row < R; row++) {
374  for (col = 0; col < row; col++)
375  ans += 2.0 * *(Aptr++) * *(Bptr++);
376  ans += *(Aptr++) * *(Bptr++); // Diagonal.
377  }
378  return ans;
379 }
380 
381 template
383 
384 template
386 
387 
388 template<typename Real>
389 Real TraceSpMat(const SpMatrix<Real> &A, const MatrixBase<Real> &B) {
390  KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols() &&
391  "KALDI_ERR: TraceSpMat: arguments have mismatched dimension");
392  MatrixIndexT R = A.NumRows();
393  Real ans = (Real)0.0;
394  const Real *Aptr = A.Data(), *Bptr = B.Data();
395  MatrixIndexT bStride = B.Stride();
396  for (MatrixIndexT r = 0;r < R;r++) {
397  for (MatrixIndexT c = 0;c < r;c++) {
398  // ans += A(r, c) * (B(r, c) + B(c, r));
399  ans += *(Aptr++) * (Bptr[r*bStride + c] + Bptr[c*bStride + r]);
400  }
401  // ans += A(r, r) * B(r, r);
402  ans += *(Aptr++) * Bptr[r*bStride + r];
403  }
404  return ans;
405 }
406 
407 template
408 float TraceSpMat(const SpMatrix<float> &A, const MatrixBase<float> &B);
409 
410 template
411 double TraceSpMat(const SpMatrix<double> &A, const MatrixBase<double> &B);
412 
413 
414 template<typename Real>
416  const SpMatrix<Real> &B, const MatrixBase<Real> &C,
417  MatrixTransposeType transC) {
418  KALDI_ASSERT((transA == kTrans?A.NumCols():A.NumRows()) ==
419  (transC == kTrans?C.NumRows():C.NumCols()) &&
420  (transA == kTrans?A.NumRows():A.NumCols()) == B.NumRows() &&
421  (transC == kTrans?C.NumCols():C.NumRows()) == B.NumRows() &&
422  "TraceMatSpMat: arguments have wrong dimension.");
423  Matrix<Real> tmp(B.NumRows(), B.NumRows());
424  tmp.AddMatMat(1.0, C, transC, A, transA, 0.0); // tmp = C * A.
425  return TraceSpMat(B, tmp);
426 }
427 
428 template
430  const SpMatrix<float> &B, const MatrixBase<float> &C,
431  MatrixTransposeType transC);
432 template
434  const SpMatrix<double> &B, const MatrixBase<double> &C,
435  MatrixTransposeType transC);
436 
437 template<typename Real>
439  const SpMatrix<Real> &B, const MatrixBase<Real> &C,
440  MatrixTransposeType transC, const SpMatrix<Real> &D) {
441  KALDI_ASSERT((transA == kTrans ?A.NumCols():A.NumRows() == D.NumCols()) &&
442  (transA == kTrans ? A.NumRows():A.NumCols() == B.NumRows()) &&
443  (transC == kTrans ? A.NumCols():A.NumRows() == B.NumCols()) &&
444  (transC == kTrans ? A.NumRows():A.NumCols() == D.NumRows()) &&
445  "KALDI_ERR: TraceMatSpMatSp: arguments have mismatched dimension.");
446  // Could perhaps optimize this more depending on dimensions of quantities.
447  Matrix<Real> tmpAB(transA == kTrans ? A.NumCols():A.NumRows(), B.NumCols());
448  tmpAB.AddMatSp(1.0, A, transA, B, 0.0);
449  Matrix<Real> tmpCD(transC == kTrans ? C.NumCols():C.NumRows(), D.NumCols());
450  tmpCD.AddMatSp(1.0, C, transC, D, 0.0);
451  return TraceMatMat(tmpAB, tmpCD, kNoTrans);
452 }
453 
454 template
456  const SpMatrix<float> &B, const MatrixBase<float> &C,
457  MatrixTransposeType transC, const SpMatrix<float> &D);
458 template
460  const SpMatrix<double> &B, const MatrixBase<double> &C,
461  MatrixTransposeType transC, const SpMatrix<double> &D);
462 
463 
464 template<typename Real>
465 bool SpMatrix<Real>::IsDiagonal(Real cutoff) const {
466  MatrixIndexT R = this->NumRows();
467  Real bad_sum = 0.0, good_sum = 0.0;
468  for (MatrixIndexT i = 0; i < R; i++) {
469  for (MatrixIndexT j = 0; j <= i; j++) {
470  if (i == j)
471  good_sum += std::abs((*this)(i, j));
472  else
473  bad_sum += std::abs((*this)(i, j));
474  }
475  }
476  return (!(bad_sum > good_sum * cutoff));
477 }
478 
479 template<typename Real>
480 bool SpMatrix<Real>::IsUnit(Real cutoff) const {
481  MatrixIndexT R = this->NumRows();
482  Real max = 0.0; // max error
483  for (MatrixIndexT i = 0; i < R; i++)
484  for (MatrixIndexT j = 0; j <= i; j++)
485  max = std::max(max, static_cast<Real>(std::abs((*this)(i, j) -
486  (i == j ? 1.0 : 0.0))));
487  return (max <= cutoff);
488 }
489 
490 template<typename Real>
491 bool SpMatrix<Real>::IsTridiagonal(Real cutoff) const {
492  MatrixIndexT R = this->NumRows();
493  Real max_abs_2diag = 0.0, max_abs_offdiag = 0.0;
494  for (MatrixIndexT i = 0; i < R; i++)
495  for (MatrixIndexT j = 0; j <= i; j++) {
496  if (j+1 < i)
497  max_abs_offdiag = std::max(max_abs_offdiag,
498  std::abs((*this)(i, j)));
499  else
500  max_abs_2diag = std::max(max_abs_2diag,
501  std::abs((*this)(i, j)));
502  }
503  return (max_abs_offdiag <= cutoff * max_abs_2diag);
504 }
505 
506 template<typename Real>
507 bool SpMatrix<Real>::IsZero(Real cutoff) const {
508  if (this->num_rows_ == 0) return true;
509  return (this->Max() <= cutoff && this->Min() >= -cutoff);
510 }
511 
512 template<typename Real>
514  Real sum = 0.0;
515  MatrixIndexT R = this->NumRows();
516  for (MatrixIndexT i = 0; i < R; i++) {
517  for (MatrixIndexT j = 0; j < i; j++)
518  sum += (*this)(i, j) * (*this)(i, j) * 2;
519  sum += (*this)(i, i) * (*this)(i, i);
520  }
521  return std::sqrt(sum);
522 }
523 
524 template<typename Real>
525 bool SpMatrix<Real>::ApproxEqual(const SpMatrix<Real> &other, float tol) const {
526  if (this->NumRows() != other.NumRows())
527  KALDI_ERR << "SpMatrix::AproxEqual, size mismatch, "
528  << this->NumRows() << " vs. " << other.NumRows();
529  SpMatrix<Real> tmp(*this);
530  tmp.AddSp(-1.0, other);
531  return (tmp.FrobeniusNorm() <= tol * std::max(this->FrobeniusNorm(), other.FrobeniusNorm()));
532 }
533 
534 // function Floor: A = Floor(B, alpha * C) ... see tutorial document.
535 template<typename Real>
537  bool verbose) {
538  MatrixIndexT dim = this->NumRows();
539  int nfloored = 0;
540  KALDI_ASSERT(C.NumRows() == dim);
541  KALDI_ASSERT(alpha > 0);
542  TpMatrix<Real> L(dim);
543  L.Cholesky(C);
544  L.Scale(std::sqrt(alpha)); // equivalent to scaling C by alpha.
545  TpMatrix<Real> LInv(L);
546  LInv.Invert();
547 
548  SpMatrix<Real> D(dim);
549  { // D = L^{-1} * (*this) * L^{-T}
550  Matrix<Real> LInvFull(LInv);
551  D.AddMat2Sp(1.0, LInvFull, kNoTrans, (*this), 0.0);
552  }
553 
554  Vector<Real> l(dim);
555  Matrix<Real> U(dim, dim);
556 
557  D.Eig(&l, &U);
558 
559  if (verbose) {
560  KALDI_LOG << "ApplyFloor: flooring following diagonal to 1: " << l;
561  }
562  for (MatrixIndexT i = 0; i < l.Dim(); i++) {
563  if (l(i) < 1.0) {
564  nfloored++;
565  l(i) = 1.0;
566  }
567  }
568  l.ApplyPow(0.5);
569  U.MulColsVec(l);
570  D.AddMat2(1.0, U, kNoTrans, 0.0);
571  { // D' := U * diag(l') * U^T ... l'=floor(l, 1)
572  Matrix<Real> LFull(L);
573  (*this).AddMat2Sp(1.0, LFull, kNoTrans, D, 0.0); // A := L * D' * L^T
574  }
575  return nfloored;
576 }
577 
578 template<typename Real>
579 Real SpMatrix<Real>::LogDet(Real *det_sign) const {
580  Real log_det;
581  SpMatrix<Real> tmp(*this);
582  // false== output not needed (saves some computation).
583  tmp.Invert(&log_det, det_sign, false);
584  return log_det;
585 }
586 
587 
588 template<typename Real>
589 int SpMatrix<Real>::ApplyFloor(Real floor) {
590  MatrixIndexT Dim = this->NumRows();
591  int nfloored = 0;
592  Vector<Real> s(Dim);
593  Matrix<Real> P(Dim, Dim);
594  (*this).Eig(&s, &P);
595  for (MatrixIndexT i = 0; i < Dim; i++) {
596  if (s(i) < floor) {
597  nfloored++;
598  s(i) = floor;
599  }
600  }
601  (*this).AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
602  return nfloored;
603 }
604 
605 template<typename Real>
606 MatrixIndexT SpMatrix<Real>::LimitCond(Real maxCond, bool invert) { // e.g. maxCond = 1.0e+05.
607  MatrixIndexT Dim = this->NumRows();
608  Vector<Real> s(Dim);
609  Matrix<Real> P(Dim, Dim);
610  (*this).SymPosSemiDefEig(&s, &P);
611  KALDI_ASSERT(maxCond > 1);
612  Real floor = s.Max() / maxCond;
613  if (floor < 0) floor = 0;
614  if (floor < 1.0e-40) {
615  KALDI_WARN << "LimitCond: limiting " << floor << " to 1.0e-40";
616  floor = 1.0e-40;
617  }
618  MatrixIndexT nfloored = 0;
619  for (MatrixIndexT i = 0; i < Dim; i++) {
620  if (s(i) <= floor) nfloored++;
621  if (invert)
622  s(i) = 1.0 / std::sqrt(std::max(s(i), floor));
623  else
624  s(i) = std::sqrt(std::max(s(i), floor));
625  }
626  P.MulColsVec(s);
627  (*this).AddMat2(1.0, P, kNoTrans, 0.0); // (*this) = P*P^T. ... (*this) = P * floor(s) * P^T ... if P was original P.
628  return nfloored;
629 }
630 
631 void SolverOptions::Check() const {
632  KALDI_ASSERT(K>10 && eps<1.0e-10);
633 }
634 
635 template<> double SolveQuadraticProblem(const SpMatrix<double> &H,
636  const VectorBase<double> &g,
637  const SolverOptions &opts,
638  VectorBase<double> *x) {
639  KALDI_ASSERT(H.NumRows() == g.Dim() && g.Dim() == x->Dim() && x->Dim() != 0);
640  opts.Check();
641  MatrixIndexT dim = x->Dim();
642  if (H.IsZero(0.0)) {
643  KALDI_WARN << "Zero quadratic term in quadratic vector problem for "
644  << opts.name << ": leaving it unchanged.";
645  return 0.0;
646  }
647  if (opts.diagonal_precondition) {
648  // We can re-cast the problem with a diagonal preconditioner to
649  // make H better-conditioned.
650  Vector<double> H_diag(dim);
651  H_diag.CopyDiagFromSp(H);
652  H_diag.ApplyFloor(std::numeric_limits<double>::min() * 1.0E+3);
653  Vector<double> H_diag_sqrt(H_diag);
654  H_diag_sqrt.ApplyPow(0.5);
655  Vector<double> H_diag_inv_sqrt(H_diag_sqrt);
656  H_diag_inv_sqrt.InvertElements();
657  Vector<double> x_scaled(*x);
658  x_scaled.MulElements(H_diag_sqrt);
659  Vector<double> g_scaled(g);
660  g_scaled.MulElements(H_diag_inv_sqrt);
661  SpMatrix<double> H_scaled(dim);
662  H_scaled.AddVec2Sp(1.0, H_diag_inv_sqrt, H, 0.0);
663  double ans;
664  SolverOptions new_opts(opts);
665  new_opts.diagonal_precondition = false;
666  ans = SolveQuadraticProblem(H_scaled, g_scaled, new_opts, &x_scaled);
667  x->CopyFromVec(x_scaled);
668  x->MulElements(H_diag_inv_sqrt);
669  return ans;
670  }
671  Vector<double> gbar(g);
672  if (opts.optimize_delta) gbar.AddSpVec(-1.0, H, *x, 1.0); // gbar = g - H x
673  Matrix<double> U(dim, dim);
674  Vector<double> l(dim);
675  H.SymPosSemiDefEig(&l, &U); // does svd H = U L V^T and checks that H == U L U^T to within a tolerance.
676  // floor l.
677  double f = std::max(static_cast<double>(opts.eps), l.Max() / opts.K);
678  MatrixIndexT nfloored = 0;
679  for (MatrixIndexT i = 0; i < dim; i++) { // floor l.
680  if (l(i) < f) {
681  nfloored++;
682  l(i) = f;
683  }
684  }
685  if (nfloored != 0 && opts.print_debug_output) {
686  KALDI_LOG << "Solving quadratic problem for " << opts.name
687  << ": floored " << nfloored<< " eigenvalues. ";
688  }
689  Vector<double> tmp(dim);
690  tmp.AddMatVec(1.0, U, kTrans, gbar, 0.0); // tmp = U^T \bar{g}
691  tmp.DivElements(l); // divide each element of tmp by l: tmp = \tilde{L}^{-1} U^T \bar{g}
692  Vector<double> delta(dim);
693  delta.AddMatVec(1.0, U, kNoTrans, tmp, 0.0); // delta = U tmp = U \tilde{L}^{-1} U^T \bar{g}
694  Vector<double> &xhat(tmp);
695  xhat.CopyFromVec(delta);
696  if (opts.optimize_delta) xhat.AddVec(1.0, *x); // xhat = x + delta.
697  double auxf_before = VecVec(g, *x) - 0.5 * VecSpVec(*x, H, *x),
698  auxf_after = VecVec(g, xhat) - 0.5 * VecSpVec(xhat, H, xhat);
699  if (auxf_after < auxf_before) { // Reject change.
700  if (auxf_after < auxf_before - 1.0e-10 && opts.print_debug_output)
701  KALDI_WARN << "Optimizing vector auxiliary function for "
702  << opts.name<< ": auxf decreased " << auxf_before
703  << " to " << auxf_after << ", change is "
704  << (auxf_after-auxf_before);
705  return 0.0;
706  } else {
707  x->CopyFromVec(xhat);
708  return auxf_after - auxf_before;
709  }
710 }
711 
712 template<> float SolveQuadraticProblem(const SpMatrix<float> &H,
713  const VectorBase<float> &g,
714  const SolverOptions &opts,
715  VectorBase<float> *x) {
716  KALDI_ASSERT(H.NumRows() == g.Dim() && g.Dim() == x->Dim() && x->Dim() != 0);
717  SpMatrix<double> Hd(H);
718  Vector<double> gd(g);
719  Vector<double> xd(*x);
720  float ans = static_cast<float>(SolveQuadraticProblem(Hd, gd, opts, &xd));
721  x->CopyFromVec(xd);
722  return ans;
723 }
724 
725 // Maximizes the auxiliary function Q(x) = tr(M^T SigmaInv Y) - 0.5 tr(SigmaInv M Q M^T).
726 // Like a numerically stable version of M := Y Q^{-1}.
727 template<typename Real>
728 Real
730  const MatrixBase<Real> &Y,
731  const SpMatrix<Real> &SigmaInv,
732  const SolverOptions &opts,
733  MatrixBase<Real> *M) {
734  KALDI_ASSERT(Q.NumRows() == M->NumCols() &&
735  SigmaInv.NumRows() == M->NumRows() && Y.NumRows() == M->NumRows()
736  && Y.NumCols() == M->NumCols() && M->NumCols() != 0);
737  opts.Check();
738  MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
739  if (Q.IsZero(0.0)) {
740  KALDI_WARN << "Zero quadratic term in quadratic matrix problem for "
741  << opts.name << ": leaving it unchanged.";
742  return 0.0;
743  }
744 
745  if (opts.diagonal_precondition) {
746  // We can re-cast the problem with a diagonal preconditioner in the space
747  // of Q (columns of M). Helps to improve the condition of Q.
748  Vector<Real> Q_diag(cols);
749  Q_diag.CopyDiagFromSp(Q);
750  Q_diag.ApplyFloor(std::numeric_limits<Real>::min() * 1.0E+3);
751  Vector<Real> Q_diag_sqrt(Q_diag);
752  Q_diag_sqrt.ApplyPow(0.5);
753  Vector<Real> Q_diag_inv_sqrt(Q_diag_sqrt);
754  Q_diag_inv_sqrt.InvertElements();
755  Matrix<Real> M_scaled(*M);
756  M_scaled.MulColsVec(Q_diag_sqrt);
757  Matrix<Real> Y_scaled(Y);
758  Y_scaled.MulColsVec(Q_diag_inv_sqrt);
759  SpMatrix<Real> Q_scaled(cols);
760  Q_scaled.AddVec2Sp(1.0, Q_diag_inv_sqrt, Q, 0.0);
761  Real ans;
762  SolverOptions new_opts(opts);
763  new_opts.diagonal_precondition = false;
764  ans = SolveQuadraticMatrixProblem(Q_scaled, Y_scaled, SigmaInv,
765  new_opts, &M_scaled);
766  M->CopyFromMat(M_scaled);
767  M->MulColsVec(Q_diag_inv_sqrt);
768  return ans;
769  }
770 
771  Matrix<Real> Ybar(Y);
772  if (opts.optimize_delta) {
773  Matrix<Real> Qfull(Q);
774  Ybar.AddMatMat(-1.0, *M, kNoTrans, Qfull, kNoTrans, 1.0);
775  } // Ybar = Y - M Q.
776  Matrix<Real> U(cols, cols);
777  Vector<Real> l(cols);
778  Q.SymPosSemiDefEig(&l, &U); // does svd Q = U L V^T and checks that Q == U L U^T to within a tolerance.
779  // floor l.
780  Real f = std::max<Real>(static_cast<Real>(opts.eps), l.Max() / opts.K);
781  MatrixIndexT nfloored = 0;
782  for (MatrixIndexT i = 0; i < cols; i++) { // floor l.
783  if (l(i) < f) { nfloored++; l(i) = f; }
784  }
785  if (nfloored != 0 && opts.print_debug_output)
786  KALDI_LOG << "Solving matrix problem for " << opts.name
787  << ": floored " << nfloored << " eigenvalues. ";
788  Matrix<Real> tmpDelta(rows, cols);
789  tmpDelta.AddMatMat(1.0, Ybar, kNoTrans, U, kNoTrans, 0.0); // tmpDelta = Ybar * U.
790  l.InvertElements(); KALDI_ASSERT(1.0/l.Max() != 0); // check not infinite. eps should take care of this.
791  tmpDelta.MulColsVec(l); // tmpDelta = Ybar * U * \tilde{L}^{-1}
792 
793  Matrix<Real> Delta(rows, cols);
794  Delta.AddMatMat(1.0, tmpDelta, kNoTrans, U, kTrans, 0.0); // Delta = Ybar * U * \tilde{L}^{-1} * U^T
795 
796  Real auxf_before, auxf_after;
797  SpMatrix<Real> MQM(rows);
798  Matrix<Real> &SigmaInvY(tmpDelta);
799  { Matrix<Real> SigmaInvFull(SigmaInv); SigmaInvY.AddMatMat(1.0, SigmaInvFull, kNoTrans, Y, kNoTrans, 0.0); }
800  { // get auxf_before. Q(x) = tr(M^T SigmaInv Y) - 0.5 tr(SigmaInv M Q M^T).
801  MQM.AddMat2Sp(1.0, *M, kNoTrans, Q, 0.0);
802  auxf_before = TraceMatMat(*M, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
803  }
804 
805  Matrix<Real> Mhat(Delta);
806  if (opts.optimize_delta) Mhat.AddMat(1.0, *M); // Mhat = Delta + M.
807 
808  { // get auxf_after.
809  MQM.AddMat2Sp(1.0, Mhat, kNoTrans, Q, 0.0);
810  auxf_after = TraceMatMat(Mhat, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
811  }
812 
813  if (auxf_after < auxf_before) {
814  if (auxf_after < auxf_before - 1.0e-10)
815  KALDI_WARN << "Optimizing matrix auxiliary function for "
816  << opts.name << ", auxf decreased "
817  << auxf_before << " to " << auxf_after << ", change is "
818  << (auxf_after-auxf_before);
819  return 0.0;
820  } else {
821  M->CopyFromMat(Mhat);
822  return auxf_after - auxf_before;
823  }
824 }
825 
826 template<typename Real>
828  const SpMatrix<Real> &P1,
829  const SpMatrix<Real> &P2,
830  const SpMatrix<Real> &Q1,
831  const SpMatrix<Real> &Q2,
832  const SolverOptions &opts,
833  MatrixBase<Real> *M) {
834  KALDI_ASSERT(Q1.NumRows() == M->NumCols() && P1.NumRows() == M->NumRows() &&
835  G.NumRows() == M->NumRows() && G.NumCols() == M->NumCols() &&
836  M->NumCols() != 0 && Q2.NumRows() == M->NumCols() &&
837  P2.NumRows() == M->NumRows());
838  MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
839  // The following check should not fail as we stipulate P1, P2 and one of Q1
840  // or Q2 must be +ve def and other Q1 or Q2 must be +ve semidef.
841  TpMatrix<Real> LInv(rows);
842  LInv.Cholesky(P1);
843  LInv.Invert(); // Will throw exception if fails.
844  SpMatrix<Real> S(rows);
845  Matrix<Real> LInvFull(LInv);
846  S.AddMat2Sp(1.0, LInvFull, kNoTrans, P2, 0.0); // S := L^{-1} P_2 L^{-T}
847  Matrix<Real> U(rows, rows);
848  Vector<Real> d(rows);
849  S.SymPosSemiDefEig(&d, &U);
850  Matrix<Real> T(rows, rows);
851  T.AddMatMat(1.0, U, kTrans, LInvFull, kNoTrans, 0.0); // T := U^T * L^{-1}
852 
853 #ifdef KALDI_PARANOID // checking mainly for errors in the code or math.
854  {
855  SpMatrix<Real> P1Trans(rows);
856  P1Trans.AddMat2Sp(1.0, T, kNoTrans, P1, 0.0);
857  KALDI_ASSERT(P1Trans.IsUnit(0.01));
858  }
859  {
860  SpMatrix<Real> P2Trans(rows);
861  P2Trans.AddMat2Sp(1.0, T, kNoTrans, P2, 0.0);
862  KALDI_ASSERT(P2Trans.IsDiagonal(0.01));
863  }
864 #endif
865 
866  Matrix<Real> TInv(T);
867  TInv.Invert();
868  Matrix<Real> Gdash(rows, cols);
869  Gdash.AddMatMat(1.0, T, kNoTrans, G, kNoTrans, 0.0); // G' = T G
870  Matrix<Real> MdashOld(rows, cols);
871  MdashOld.AddMatMat(1.0, TInv, kTrans, *M, kNoTrans, 0.0); // M' = T^{-T} M
872  Matrix<Real> MdashNew(MdashOld);
873  Real objf_impr = 0.0;
874  for (MatrixIndexT n = 0; n < rows; n++) {
875  SpMatrix<Real> Qsum(Q1);
876  Qsum.AddSp(d(n), Q2);
877  SubVector<Real> mdash_n = MdashNew.Row(n);
878  SubVector<Real> gdash_n = Gdash.Row(n);
879 
880  Matrix<Real> QsumInv(Qsum);
881  try {
882  QsumInv.Invert();
883  Real old_objf = VecVec(mdash_n, gdash_n)
884  - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
885  mdash_n.AddMatVec(1.0, QsumInv, kNoTrans, gdash_n, 0.0); // m'_n := g'_n * (Q_1 + d_n Q_2)^{-1}
886  Real new_objf = VecVec(mdash_n, gdash_n)
887  - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
888  if (new_objf < old_objf) {
889  if (new_objf < old_objf - 1.0e-05) {
890  KALDI_WARN << "In double quadratic matrix problem: objective "
891  "function decreasing during optimization of " << opts.name
892  << ", " << old_objf << "->" << new_objf << ", change is "
893  << (new_objf - old_objf);
894  KALDI_ERR << "Auxiliary function decreasing."; // Will be caught.
895  } else { // Reset to old value, didn't improve (very close to optimum).
896  MdashNew.Row(n).CopyFromVec(MdashOld.Row(n));
897  }
898  }
899  objf_impr += new_objf - old_objf;
900  }
901  catch (...) {
902  KALDI_WARN << "Matrix inversion or optimization failed during double "
903  "quadratic problem, solving for" << opts.name
904  << ": trying more stable approach.";
905  objf_impr += SolveQuadraticProblem(Qsum, gdash_n, opts, &mdash_n);
906  }
907  }
908  M->AddMatMat(1.0, T, kTrans, MdashNew, kNoTrans, 0.0); // M := T^T M'.
909  return objf_impr;
910 }
911 
912 // rank-one update, this <-- this + alpha V V'
913 template<>
914 template<>
915 void SpMatrix<float>::AddVec2(const float alpha, const VectorBase<float> &v) {
916  KALDI_ASSERT(v.Dim() == this->NumRows());
917  cblas_Xspr(v.Dim(), alpha, v.Data(), 1,
918  this->data_);
919 }
920 
921 template<class Real>
922 void SpMatrix<Real>::AddVec2Sp(const Real alpha, const VectorBase<Real> &v,
923  const SpMatrix<Real> &S, const Real beta) {
924  KALDI_ASSERT(v.Dim() == this->NumRows() && S.NumRows() == this->NumRows());
925  const Real *Sdata = S.Data();
926  const Real *vdata = v.Data();
927  Real *data = this->data_;
928  MatrixIndexT dim = this->num_rows_;
929  for (MatrixIndexT r = 0; r < dim; r++)
930  for (MatrixIndexT c = 0; c <= r; c++, Sdata++, data++)
931  *data = beta * *data + alpha * vdata[r] * vdata[c] * *Sdata;
932 }
933 
934 
935 // rank-one update, this <-- this + alpha V V'
936 template<>
937 template<>
938 void SpMatrix<double>::AddVec2(const double alpha, const VectorBase<double> &v) {
939  KALDI_ASSERT(v.Dim() == num_rows_);
940  cblas_Xspr(v.Dim(), alpha, v.Data(), 1, data_);
941 }
942 
943 
944 template<typename Real>
945 template<typename OtherReal>
946 void SpMatrix<Real>::AddVec2(const Real alpha, const VectorBase<OtherReal> &v) {
947  KALDI_ASSERT(v.Dim() == this->NumRows());
948  Real *data = this->data_;
949  const OtherReal *v_data = v.Data();
950  MatrixIndexT nr = this->num_rows_;
951  for (MatrixIndexT i = 0; i < nr; i++)
952  for (MatrixIndexT j = 0; j <= i; j++, data++)
953  *data += alpha * v_data[i] * v_data[j];
954 }
955 
956 // instantiate the template above.
957 template
958 void SpMatrix<float>::AddVec2(const float alpha, const VectorBase<double> &v);
959 template
960 void SpMatrix<double>::AddVec2(const double alpha, const VectorBase<float> &v);
961 
962 
963 template<typename Real>
964 Real VecSpVec(const VectorBase<Real> &v1, const SpMatrix<Real> &M,
965  const VectorBase<Real> &v2) {
966  MatrixIndexT D = M.NumRows();
967  KALDI_ASSERT(v1.Dim() == D && v1.Dim() == v2.Dim());
968  Vector<Real> tmp_vec(D);
969  cblas_Xspmv(D, 1.0, M.Data(), v1.Data(), 1, 0.0, tmp_vec.Data(), 1);
970  return VecVec(tmp_vec, v2);
971 }
972 
973 template
974 float VecSpVec(const VectorBase<float> &v1, const SpMatrix<float> &M,
975  const VectorBase<float> &v2);
976 template
977 double VecSpVec(const VectorBase<double> &v1, const SpMatrix<double> &M,
978  const VectorBase<double> &v2);
979 
980 
981 template<typename Real>
983  const Real alpha, const MatrixBase<Real> &M,
984  MatrixTransposeType transM, const SpMatrix<Real> &A, const Real beta) {
985  if (transM == kNoTrans) {
986  KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
987  } else {
988  KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
989  }
990  Vector<Real> tmp_vec(A.NumRows());
991  Real *tmp_vec_data = tmp_vec.Data();
992  SpMatrix<Real> tmp_A;
993  const Real *p_A_data = A.Data();
994  Real *p_row_data = this->Data();
995  MatrixIndexT M_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows()),
996  M_same_dim = (transM == kNoTrans ? M.NumRows() : M.NumCols()),
997  M_stride = M.Stride(), dim = this->NumRows();
998  KALDI_ASSERT(M_same_dim == dim);
999 
1000  const Real *M_data = M.Data();
1001 
1002  if (this->Data() <= A.Data() + A.SizeInBytes() &&
1003  this->Data() + this->SizeInBytes() >= A.Data()) {
1004  // Matrices A and *this overlap. Make copy of A
1005  tmp_A.Resize(A.NumRows());
1006  tmp_A.CopyFromSp(A);
1007  p_A_data = tmp_A.Data();
1008  }
1009 
1010  if (transM == kNoTrans) {
1011  for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
1012  cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.RowData(r), 1, 0.0, tmp_vec_data, 1);
1013  cblas_Xgemv(transM, r+1, M_other_dim, alpha, M_data, M_stride,
1014  tmp_vec_data, 1, beta, p_row_data, 1);
1015  }
1016  } else {
1017  for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
1018  cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.Data() + r, M.Stride(), 0.0, tmp_vec_data, 1);
1019  cblas_Xgemv(transM, M_other_dim, r+1, alpha, M_data, M_stride,
1020  tmp_vec_data, 1, beta, p_row_data, 1);
1021  }
1022  }
1023 }
1024 
1025 template<typename Real>
1027  const Real alpha, const MatrixBase<Real> &M,
1028  MatrixTransposeType transM, const SpMatrix<Real> &A,
1029  const Real beta) {
1030  KALDI_ASSERT((transM == kNoTrans && M.NumCols() == A.NumRows()) ||
1031  (transM == kTrans && M.NumRows() == A.NumRows()));
1032  if (transM == kNoTrans) {
1033  KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
1034  } else {
1035  KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
1036  }
1037  MatrixIndexT Adim = A.NumRows(), dim = this->num_rows_;
1038 
1039  Matrix<Real> temp_A(A); // represent A as full matrix.
1040  Matrix<Real> temp_MA(dim, Adim);
1041  temp_MA.AddSmatMat(1.0, M, transM, temp_A, kNoTrans, 0.0);
1042 
1043  // Next-- we want to do *this = alpha * temp_MA * M^T + beta * *this.
1044  // To make it sparse vector multiplies, since M is sparse, we'd like
1045  // to do: for each column c, (*this column c) += temp_MA * (M^T's column c.)
1046  // [ignoring the alpha and beta here.]
1047  // It's not convenient to process columns in the symmetric
1048  // packed format because they don't have a constant stride. However,
1049  // we can use the fact that temp_MA * M is symmetric, to just assign
1050  // each row of *this instead of each column.
1051  // So the final iteration is:
1052  // for i = 0... dim-1,
1053  // [the i'th row of *this] = beta * [the i'th row of *this] + alpha *
1054  // temp_MA * [the i'th column of M].
1055  // Of course, we only process the first 0 ... i elements of this row,
1056  // as that's all that are kept in the symmetric packed format.
1057 
1058  Matrix<Real> temp_this(*this);
1059  Real *data = this->data_;
1060  const Real *Mdata = M.Data(), *MAdata = temp_MA.Data();
1061  MatrixIndexT temp_MA_stride = temp_MA.Stride(), Mstride = M.Stride();
1062 
1063  if (transM == kNoTrans) {
1064  // The column of M^T corresponds to the rows of the supplied matrix.
1065  for (MatrixIndexT i = 0; i < dim; i++, data += i) {
1066  MatrixIndexT num_rows = i + 1, num_cols = Adim;
1067  Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
1068  temp_MA_stride, Mdata + (i * Mstride), 1, beta, data, 1);
1069  }
1070  } else {
1071  // The column of M^T corresponds to the columns of the supplied matrix.
1072  for (MatrixIndexT i = 0; i < dim; i++, data += i) {
1073  MatrixIndexT num_rows = i + 1, num_cols = Adim;
1074  Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
1075  temp_MA_stride, Mdata + i, Mstride, beta, data, 1);
1076  }
1077  }
1078 }
1079 
1080 template<typename Real>
1081 void SpMatrix<Real>::AddMat2Vec(const Real alpha,
1082  const MatrixBase<Real> &M,
1083  MatrixTransposeType transM,
1084  const VectorBase<Real> &v,
1085  const Real beta) {
1086  this->Scale(beta);
1087  KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows() &&
1088  M.NumCols() == v.Dim()) ||
1089  (transM == kTrans && this->NumRows() == M.NumCols() &&
1090  M.NumRows() == v.Dim()));
1091 
1092  if (transM == kNoTrans) {
1093  const Real *Mdata = M.Data(), *vdata = v.Data();
1094  Real *data = this->data_;
1095  MatrixIndexT dim = this->NumRows(), mcols = M.NumCols(),
1096  mstride = M.Stride();
1097  for (MatrixIndexT col = 0; col < mcols; col++, vdata++, Mdata += 1)
1098  cblas_Xspr(dim, *vdata*alpha, Mdata, mstride, data);
1099  } else {
1100  const Real *Mdata = M.Data(), *vdata = v.Data();
1101  Real *data = this->data_;
1102  MatrixIndexT dim = this->NumRows(), mrows = M.NumRows(),
1103  mstride = M.Stride();
1104  for (MatrixIndexT row = 0; row < mrows; row++, vdata++, Mdata += mstride)
1105  cblas_Xspr(dim, *vdata*alpha, Mdata, 1, data);
1106  }
1107 }
1108 
1109 template<typename Real>
1110 void SpMatrix<Real>::AddMat2(const Real alpha, const MatrixBase<Real> &M,
1111  MatrixTransposeType transM, const Real beta) {
1112  KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows())
1113  || (transM == kTrans && this->NumRows() == M.NumCols()));
1114 
1115  // Cblas has no function *sprk (i.e. symmetric packed rank-k update), so we
1116  // use as temporary storage a regular matrix of which we only access its lower
1117  // triangle
1118 
1119  MatrixIndexT this_dim = this->NumRows(),
1120  m_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows());
1121 
1122  if (this_dim == 0) return;
1123  if (alpha == 0.0) {
1124  if (beta != 1.0) this->Scale(beta);
1125  return;
1126  }
1127 
1128  Matrix<Real> temp_mat(*this); // wastefully copies upper triangle too, but this
1129  // doesn't dominate O(N) time.
1130 
1131  // This function call is hard-coded to update the lower triangle.
1132  cblas_Xsyrk(transM, this_dim, m_other_dim, alpha, M.Data(),
1133  M.Stride(), beta, temp_mat.Data(), temp_mat.Stride());
1134 
1135  this->CopyFromMat(temp_mat, kTakeLower);
1136 }
1137 
1138 template<typename Real>
1139 void SpMatrix<Real>::AddTp2Sp(const Real alpha, const TpMatrix<Real> &T,
1140  MatrixTransposeType transM, const SpMatrix<Real> &A,
1141  const Real beta) {
1142  Matrix<Real> Tmat(T);
1143  AddMat2Sp(alpha, Tmat, transM, A, beta);
1144 }
1145 
1146 template<typename Real>
1147 void SpMatrix<Real>::AddVecVec(const Real alpha, const VectorBase<Real> &v,
1148  const VectorBase<Real> &w) {
1149  int32 dim = this->NumRows();
1150  KALDI_ASSERT(dim == v.Dim() && dim == w.Dim() && dim > 0);
1151  cblas_Xspr2(dim, alpha, v.Data(), 1, w.Data(), 1, this->data_);
1152 }
1153 
1154 
1155 template<typename Real>
1156 void SpMatrix<Real>::AddTp2(const Real alpha, const TpMatrix<Real> &T,
1157  MatrixTransposeType transM, const Real beta) {
1158  Matrix<Real> Tmat(T);
1159  AddMat2(alpha, Tmat, transM, beta);
1160 }
1161 
1162 
1163 // Explicit instantiation of the class.
1164 // This needs to be after the definition of all the class member functions.
1165 
1166 template class SpMatrix<float>;
1167 template class SpMatrix<double>;
1168 
1169 
1170 template<typename Real>
1172  MatrixIndexT adim = A.NumRows();
1173  KALDI_ASSERT(adim == B.NumRows());
1174  MatrixIndexT dim = (adim*(adim+1))/2;
1175  return cblas_Xdot(dim, A.Data(), 1, B.Data(), 1);
1176 }
1177 // Instantiate the template above.
1178 template
1179 double TraceSpSpLower(const SpMatrix<double> &A, const SpMatrix<double> &B);
1180 template
1181 float TraceSpSpLower(const SpMatrix<float> &A, const SpMatrix<float> &B);
1182 
1183 // Instantiate the template above.
1184 template float SolveQuadraticMatrixProblem(const SpMatrix<float> &Q,
1185  const MatrixBase<float> &Y,
1186  const SpMatrix<float> &SigmaInv,
1187  const SolverOptions &opts,
1188  MatrixBase<float> *M);
1189 template double SolveQuadraticMatrixProblem(const SpMatrix<double> &Q,
1190  const MatrixBase<double> &Y,
1191  const SpMatrix<double> &SigmaInv,
1192  const SolverOptions &opts,
1193  MatrixBase<double> *M);
1194 
1195 // Instantiate the template above.
1196 template float SolveDoubleQuadraticMatrixProblem(
1197  const MatrixBase<float> &G,
1198  const SpMatrix<float> &P1,
1199  const SpMatrix<float> &P2,
1200  const SpMatrix<float> &Q1,
1201  const SpMatrix<float> &Q2,
1202  const SolverOptions &opts,
1203  MatrixBase<float> *M);
1204 
1205 template double SolveDoubleQuadraticMatrixProblem(
1206  const MatrixBase<double> &G,
1207  const SpMatrix<double> &P1,
1208  const SpMatrix<double> &P2,
1209  const SpMatrix<double> &Q1,
1210  const SpMatrix<double> &Q2,
1211  const SolverOptions &opts,
1212  MatrixBase<double> *M);
1213 
1214 
1215 
1216 } // namespace kaldi
void AddMat2(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const Real beta)
rank-N update: if (transM == kNoTrans) (*this) = beta*(*this) + alpha * M * M^T, or (if transM == kTr...
Definition: sp-matrix.cc:1110
This code computes Goodness of Pronunciation (GOP) and extracts phone-level pronunciation feature for...
Definition: chain.dox:20
void AddVec2Sp(const Real alpha, const VectorBase< Real > &v, const SpMatrix< Real > &S, const Real beta)
Does *this = beta * *thi + alpha * diag(v) * S * diag(v)
Definition: sp-matrix.cc:922
bool IsUnit(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:480
This class provides a way for switching between double and float types.
Definition: matrix-common.h:84
bool IsDiagonal(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:465
Packed symetric matrix class.
Definition: matrix-common.h:62
void cblas_Xsyrk(const MatrixTransposeType trans, const MatrixIndexT dim_c, const MatrixIndexT other_dim_a, const float alpha, const float *A, const MatrixIndexT a_stride, const float beta, float *C, const MatrixIndexT c_stride)
This class describes the options for maximizing various quadratic objective functions.
Definition: sp-matrix.h:443
void Scale(Real c)
double SolveQuadraticProblem(const SpMatrix< double > &H, const VectorBase< double > &g, const SolverOptions &opts, VectorBase< double > *x)
Definition: sp-matrix.cc:635
void clapack_Xsptri(KaldiBlasInt *num_rows, float *Mdata, KaldiBlasInt *ipiv, float *work, KaldiBlasInt *result)
bool ApproxEqual(const SpMatrix< Real > &other, float tol=0.01) const
Returns true if ((*this)-other).FrobeniusNorm() <= tol*(*this).FrobeniusNorma()
Definition: sp-matrix.cc:525
void cblas_Xspr2(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, const float *Ydata, MatrixIndexT incY, float *Adata)
void AddVecVec(const Real alpha, const VectorBase< Real > &v, const VectorBase< Real > &w)
rank-two update, this <– this + alpha (v w&#39; + w v&#39;).
Definition: sp-matrix.cc:1147
void AddMat2Vec(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const VectorBase< Real > &v, const Real beta=0.0)
Extension of rank-N update: this <– beta*this + alpha * M * diag(v) * M^T.
Definition: sp-matrix.cc:1081
MatrixIndexT NumCols() const
Returns number of columns (or zero for empty matrix).
Definition: kaldi-matrix.h:67
Real SolveQuadraticMatrixProblem(const SpMatrix< Real > &Q, const MatrixBase< Real > &Y, const SpMatrix< Real > &SigmaInv, const SolverOptions &opts, MatrixBase< Real > *M)
Maximizes the auxiliary function : Like a numerically stable version of .
Definition: sp-matrix.cc:729
Base class which provides matrix operations not involving resizing or allocation. ...
Definition: kaldi-matrix.h:49
const Real * Data() const
Gives pointer to raw data (const).
Definition: kaldi-matrix.h:79
Real Trace() const
Definition: sp-matrix.cc:171
void InvertDouble(Real *logdet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
Definition: sp-matrix.cc:312
void Xgemv_sparsevec(MatrixTransposeType trans, MatrixIndexT num_rows, MatrixIndexT num_cols, Real alpha, const Real *Mdata, MatrixIndexT stride, const Real *xdata, MatrixIndexT incX, Real beta, Real *ydata, MatrixIndexT incY)
Real FrobeniusNorm() const
sqrt of sum of square elements.
Definition: sp-matrix.cc:513
Real TraceSpMat(const SpMatrix< Real > &A, const MatrixBase< Real > &B)
Returns tr(A B).
Definition: sp-matrix.cc:389
Real * RowData(MatrixIndexT i)
Returns pointer to data for one row (non-const)
Definition: kaldi-matrix.h:87
void AddMat(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transA=kNoTrans)
*this += alpha * M [or M^T]
void swap(basic_filebuf< CharT, Traits > &x, basic_filebuf< CharT, Traits > &y)
kaldi::int32 int32
void ApplyPow(Real exponent)
Takes matrix to a fraction power via Svd.
Definition: sp-matrix.cc:91
bool IsPosDef() const
Definition: sp-matrix.cc:75
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
A class for storing matrices.
Definition: kaldi-matrix.h:823
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.
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.
bool IsTridiagonal(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:491
MatrixIndexT NumRows() const
uint64 data_
Real SolveDoubleQuadraticMatrixProblem(const MatrixBase< Real > &G, const SpMatrix< Real > &P1, const SpMatrix< Real > &P2, const SpMatrix< Real > &Q1, const SpMatrix< Real > &Q2, const SolverOptions &opts, MatrixBase< Real > *M)
Maximizes the auxiliary function : Encountered in matrix update with a prior.
Definition: sp-matrix.cc:827
void CopyFromSp(const SpMatrix< Real > &other)
Definition: sp-matrix.h:85
MatrixIndexT num_rows_
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
void CopyFromSp(const SpMatrix< OtherReal > &M)
Copy given spmatrix. (no resize is done).
int ApplyFloor(const SpMatrix< Real > &Floor, Real alpha=1.0, bool verbose=false)
Floors this symmetric matrix to the matrix alpha * Floor, where the matrix Floor is positive definite...
Definition: sp-matrix.cc:536
void ApplyFloor(Real floor_val, MatrixIndexT *floored_count=nullptr)
Applies floor to all elements.
Definition: kaldi-vector.h:149
float cblas_Xdot(const int N, const float *const X, const int incX, const float *const Y, const int incY)
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
MatrixIndexT NumCols() const
void Cholesky(const SpMatrix< Real > &orig)
Definition: tp-matrix.cc:88
MatrixIndexT Stride() const
Stride (distance in memory between each row). Will be >= NumCols.
Definition: kaldi-matrix.h:70
int32 MatrixIndexT
Definition: matrix-common.h:98
void AddSmat2Sp(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const SpMatrix< Real > &A, const Real beta=0.0)
This is a version of AddMat2Sp specialized for when M is fairly sparse.
Definition: sp-matrix.cc:1026
Real LogPosDefDet() const
Computes log determinant but only for +ve-def matrices (it uses Cholesky).
Definition: sp-matrix.cc:36
Real TraceMatSpMat(const MatrixBase< Real > &A, MatrixTransposeType transA, const SpMatrix< Real > &B, const MatrixBase< Real > &C, MatrixTransposeType transC)
Returns tr(A B C) (A and C may be transposed as specified by transA and transC).
Definition: sp-matrix.cc:415
const SubVector< Real > Row(MatrixIndexT i) const
Return specific row of matrix [const].
Definition: kaldi-matrix.h:188
double Log(double x)
Definition: kaldi-math.h:100
void Check() const
Definition: sp-matrix.cc:631
MatrixIndexT LimitCond(Real maxCond=1.0e+5, bool invert=false)
Definition: sp-matrix.cc:606
void MulElements(const VectorBase< Real > &v)
Multiply element-by-element by another vector.
void AddSp(const Real alpha, const SpMatrix< Real > &Ma)
Definition: sp-matrix.h:211
#define KALDI_MEMALIGN(align, size, pp_orig)
Definition: kaldi-utils.h:58
struct rnnlm::@11::@12 n
void AddTp2(const Real alpha, const TpMatrix< Real > &T, MatrixTransposeType transM, const Real beta=0.0)
The following function does: this <– beta*this + alpha * T * T^T.
Definition: sp-matrix.cc:1156
void AddMatSp(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const SpMatrix< Real > &B, const Real beta)
this <– beta*this + alpha*A*B.
Definition: kaldi-matrix.h:708
void AddMatMat(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const Real beta)
template float TraceSpSp< float, double >(const SpMatrix< float > &A, const SpMatrix< double > &B)
std::string name
Definition: sp-matrix.h:446
#define KALDI_ERR
Definition: kaldi-error.h:147
Real Max() const
Returns the maximum value of any element, or -infinity for the empty vector.
double TraceSpSp(const SpMatrix< double > &A, const SpMatrix< double > &B)
Definition: sp-matrix.cc:326
#define KALDI_MEMALIGN_FREE(x)
Definition: kaldi-utils.h:60
void cblas_Xspr(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, float *Adata)
Real VecSpVec(const VectorBase< Real > &v1, const SpMatrix< Real > &M, const VectorBase< Real > &v2)
Computes v1^T * M * v2.
Definition: sp-matrix.cc:964
Real TraceSpSpLower(const SpMatrix< Real > &A, const SpMatrix< Real > &B)
Definition: sp-matrix.cc:1171
Packed symetric matrix class.
Definition: matrix-common.h:63
#define KALDI_WARN
Definition: kaldi-error.h:150
Real TraceMatMat(const MatrixBase< Real > &A, const MatrixBase< Real > &B, MatrixTransposeType trans)
We need to declare this here as it will be a friend function.
Real * Data()
Returns a pointer to the start of the vector&#39;s data.
Definition: kaldi-vector.h:70
MatrixIndexT Dim() const
Returns the dimension of the vector.
Definition: kaldi-vector.h:64
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
Real TraceMatSpMatSp(const MatrixBase< Real > &A, MatrixTransposeType transA, const SpMatrix< Real > &B, const MatrixBase< Real > &C, MatrixTransposeType transC, const SpMatrix< Real > &D)
Returns tr (A B C D) (A and C may be transposed as specified by transA and transB).
Definition: sp-matrix.cc:438
void Swap(SpMatrix *other)
Shallow swap.
Definition: sp-matrix.cc:51
void MulColsVec(const VectorBase< Real > &scale)
Equivalent to (*this) = (*this) * diag(scale).
void AddSmatMat(Real alpha, const SparseMatrix< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, Real beta)
(*this) = alpha * op(A) * B + beta * (*this), where A is sparse.
void CopyFromMat(const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
Definition: sp-matrix.cc:112
template double TraceSpSp< double, float >(const SpMatrix< double > &A, const SpMatrix< float > &B)
A class representing a vector.
Definition: kaldi-vector.h:406
void InvertElements()
Invert all elements.
size_t SizeInBytes() const
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
MatrixIndexT NumRows() const
Returns number of rows (or zero for empty matrix).
Definition: kaldi-matrix.h:64
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
void cblas_Xgemv(MatrixTransposeType trans, MatrixIndexT num_rows, MatrixIndexT num_cols, float alpha, const float *Mdata, MatrixIndexT stride, const float *xdata, MatrixIndexT incX, float beta, float *ydata, MatrixIndexT incY)
MatrixTransposeType
Definition: matrix-common.h:32
void DivElements(const VectorBase< Real > &v)
Divide element-by-element by a vector.
void cblas_Xspmv(const float alpha, const int num_rows, const float *Mdata, const float *v, const int v_inc, const float beta, float *y, const int y_inc)
void Resize(MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
Definition: sp-matrix.h:81
void Invert(Real *log_det=NULL, Real *det_sign=NULL, bool inverse_needed=true)
matrix inverse.
Definition: kaldi-matrix.cc:38
void SymPosSemiDefEig(VectorBase< Real > *s, MatrixBase< Real > *P, Real tolerance=0.001) const
This is the version of SVD that we implement for symmetric positive definite matrices.
Definition: sp-matrix.cc:57
Provides a vector abstraction class.
Definition: kaldi-vector.h:41
void clapack_Xsptrf(KaldiBlasInt *num_rows, float *Mdata, KaldiBlasInt *ipiv, KaldiBlasInt *result)
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) ...
bool IsZero(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:507
Real MaxAbsEig() const
Returns the maximum of the absolute values of any of the eigenvalues.
Definition: sp-matrix.cc:67
void CopyDiagFromSp(const SpMatrix< Real > &M)
Extracts the diagonal of a symmetric matrix.
Definition: kaldi-vector.h:309
void AddDiagVec(const Real alpha, const VectorBase< OtherReal > &v)
diagonal update, this <– this + diag(v)
Definition: sp-matrix.cc:183
Represents a non-allocating general vector which can be defined as a sub-vector of higher-level vecto...
Definition: kaldi-vector.h:501
Real LogDet(Real *det_sign=NULL) const
Definition: sp-matrix.cc:579