sp-matrix.h
Go to the documentation of this file.
1 // matrix/sp-matrix.h
2 
3 // Copyright 2009-2011 Ondrej Glembek; Microsoft Corporation; Lukas Burget;
4 // Saarland University; Ariya Rastrow; Yanmin Qian;
5 // Jan Silovsky
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 #ifndef KALDI_MATRIX_SP_MATRIX_H_
22 #define KALDI_MATRIX_SP_MATRIX_H_
23 
24 #include <algorithm>
25 #include <vector>
26 
27 #include "matrix/packed-matrix.h"
28 
29 namespace kaldi {
30 
31 
34 template<typename Real> class SpMatrix;
35 
36 
40 template<typename Real>
41 class SpMatrix : public PackedMatrix<Real> {
42  friend class CuSpMatrix<Real>;
43  public:
44  // so it can use our assignment operator.
45  friend class std::vector<Matrix<Real> >;
46 
47  SpMatrix(): PackedMatrix<Real>() {}
48 
51 
52  explicit SpMatrix(const CuSpMatrix<Real> &cu);
53 
54  explicit SpMatrix(MatrixIndexT r, MatrixResizeType resize_type = kSetZero)
55  : PackedMatrix<Real>(r, resize_type) {}
56 
57  SpMatrix(const SpMatrix<Real> &orig)
58  : PackedMatrix<Real>(orig) {}
59 
60  template<typename OtherReal>
61  explicit SpMatrix(const SpMatrix<OtherReal> &orig)
62  : PackedMatrix<Real>(orig) {}
63 
64 #ifdef KALDI_PARANOID
65  explicit SpMatrix(const MatrixBase<Real> & orig,
66  SpCopyType copy_type = kTakeMeanAndCheck)
68  CopyFromMat(orig, copy_type);
69  }
70 #else
71  explicit SpMatrix(const MatrixBase<Real> & orig,
72  SpCopyType copy_type = kTakeMean)
73  : PackedMatrix<Real>(orig.NumRows(), kUndefined) {
74  CopyFromMat(orig, copy_type);
75  }
76 #endif
77 
79  void Swap(SpMatrix *other);
80 
81  inline void Resize(MatrixIndexT nRows, MatrixResizeType resize_type = kSetZero) {
82  PackedMatrix<Real>::Resize(nRows, resize_type);
83  }
84 
85  void CopyFromSp(const SpMatrix<Real> &other) {
87  }
88 
89  template<typename OtherReal>
90  void CopyFromSp(const SpMatrix<OtherReal> &other) {
92  }
93 
94 #ifdef KALDI_PARANOID
95  void CopyFromMat(const MatrixBase<Real> &orig,
96  SpCopyType copy_type = kTakeMeanAndCheck);
97 #else // different default arg if non-paranoid mode.
98  void CopyFromMat(const MatrixBase<Real> &orig,
99  SpCopyType copy_type = kTakeMean);
100 #endif
101 
102  inline Real operator() (MatrixIndexT r, MatrixIndexT c) const {
103  // if column is less than row, then swap these as matrix is stored
104  // as upper-triangular... only allowed for const matrix object.
105  if (static_cast<UnsignedMatrixIndexT>(c) >
106  static_cast<UnsignedMatrixIndexT>(r))
107  std::swap(c, r);
108  // c<=r now so don't have to check c.
109  KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
110  static_cast<UnsignedMatrixIndexT>(this->num_rows_));
111  return *(this->data_ + (r*(r+1)) / 2 + c);
112  // Duplicating code from PackedMatrix.h
113  }
114 
115  inline Real &operator() (MatrixIndexT r, MatrixIndexT c) {
116  if (static_cast<UnsignedMatrixIndexT>(c) >
117  static_cast<UnsignedMatrixIndexT>(r))
118  std::swap(c, r);
119  // c<=r now so don't have to check c.
120  KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
121  static_cast<UnsignedMatrixIndexT>(this->num_rows_));
122  return *(this->data_ + (r * (r + 1)) / 2 + c);
123  // Duplicating code from PackedMatrix.h
124  }
125 
128  return *this;
129  }
130 
132 
136  void Invert(Real *logdet = NULL, Real *det_sign= NULL,
137  bool inverse_needed = true);
138 
139  // Below routine does inversion in double precision,
140  // even for single-precision object.
141  void InvertDouble(Real *logdet = NULL, Real *det_sign = NULL,
142  bool inverse_needed = true);
143 
145  inline Real Cond() const {
146  Matrix<Real> tmp(*this);
147  return tmp.Cond();
148  }
149 
153  void ApplyPow(Real exponent);
154 
161  Real tolerance = 0.001) const;
162 
169  void Eig(VectorBase<Real> *s, MatrixBase<Real> *P = NULL) const;
170 
196  MatrixIndexT lanczos_dim = 0) const;
197 
198 
201  Real MaxAbsEig() const;
202 
203  void PrintEigs(const char *name) {
204  Vector<Real> s((*this).NumRows());
205  Matrix<Real> P((*this).NumRows(), (*this).NumCols());
206  SymPosSemiDefEig(&s, &P);
207  KALDI_LOG << "PrintEigs: " << name << ": " << s;
208  }
209 
210  bool IsPosDef() const; // returns true if Cholesky succeeds.
211  void AddSp(const Real alpha, const SpMatrix<Real> &Ma) {
212  this->AddPacked(alpha, Ma);
213  }
214 
219  Real LogPosDefDet() const;
220 
221  Real LogDet(Real *det_sign = NULL) const;
222 
224  template<typename OtherReal>
225  void AddVec2(const Real alpha, const VectorBase<OtherReal> &v);
226 
228  void AddVecVec(const Real alpha, const VectorBase<Real> &v,
229  const VectorBase<Real> &w);
230 
232  void AddVec2Sp(const Real alpha, const VectorBase<Real> &v,
233  const SpMatrix<Real> &S, const Real beta);
234 
236  template<typename OtherReal>
237  void AddDiagVec(const Real alpha, const VectorBase<OtherReal> &v);
238 
245  void AddMat2(const Real alpha, const MatrixBase<Real> &M,
246  MatrixTransposeType transM, const Real beta);
247 
252  void AddMat2Sp(const Real alpha, const MatrixBase<Real> &M,
253  MatrixTransposeType transM, const SpMatrix<Real> &A,
254  const Real beta = 0.0);
255 
258  void AddSmat2Sp(const Real alpha, const MatrixBase<Real> &M,
259  MatrixTransposeType transM, const SpMatrix<Real> &A,
260  const Real beta = 0.0);
261 
268  void AddTp2Sp(const Real alpha, const TpMatrix<Real> &T,
269  MatrixTransposeType transM, const SpMatrix<Real> &A,
270  const Real beta = 0.0);
271 
278  void AddTp2(const Real alpha, const TpMatrix<Real> &T,
279  MatrixTransposeType transM, const Real beta = 0.0);
280 
285  void AddMat2Vec(const Real alpha, const MatrixBase<Real> &M,
286  MatrixTransposeType transM, const VectorBase<Real> &v,
287  const Real beta = 0.0);
288 
289 
298  int ApplyFloor(const SpMatrix<Real> &Floor, Real alpha = 1.0,
299  bool verbose = false);
300 
306  int ApplyFloor(Real floor);
307 
308  bool IsDiagonal(Real cutoff = 1.0e-05) const;
309  bool IsUnit(Real cutoff = 1.0e-05) const;
310  bool IsZero(Real cutoff = 1.0e-05) const;
311  bool IsTridiagonal(Real cutoff = 1.0e-05) const;
312 
314  Real FrobeniusNorm() const;
315 
318  bool ApproxEqual(const SpMatrix<Real> &other, float tol = 0.01) const;
319 
320  // LimitCond:
321  // Limits the condition of symmetric positive semidefinite matrix to
322  // a specified value
323  // by flooring all eigenvalues to a positive number which is some multiple
324  // of the largest one (or zero if there are no positive eigenvalues).
325  // Takes the condition number we are willing to accept, and floors
326  // eigenvalues to the largest eigenvalue divided by this.
327  // Returns #eigs floored or already equal to the floor.
328  // Throws exception if input is not positive definite.
329  // returns #floored.
330  MatrixIndexT LimitCond(Real maxCond = 1.0e+5, bool invert = false);
331 
332  // as LimitCond but all done in double precision. // returns #floored.
333  MatrixIndexT LimitCondDouble(Real maxCond = 1.0e+5, bool invert = false) {
334  SpMatrix<double> dmat(*this);
335  MatrixIndexT ans = dmat.LimitCond(maxCond, invert);
336  (*this).CopyFromSp(dmat);
337  return ans;
338  }
339  Real Trace() const;
340 
346 
356  void Qr(MatrixBase<Real> *Q);
357 
358  private:
360  Real tolerance, int recurse) const;
361 };
362 
364 
367 
368 
370 float TraceSpSp(const SpMatrix<float> &A, const SpMatrix<float> &B);
371 double TraceSpSp(const SpMatrix<double> &A, const SpMatrix<double> &B);
372 
373 
374 template<typename Real>
375 inline bool ApproxEqual(const SpMatrix<Real> &A,
376  const SpMatrix<Real> &B, Real tol = 0.01) {
377  return A.ApproxEqual(B, tol);
378 }
379 
380 template<typename Real>
381 inline void AssertEqual(const SpMatrix<Real> &A,
382  const SpMatrix<Real> &B, Real tol = 0.01) {
383  KALDI_ASSERT(ApproxEqual(A, B, tol));
384 }
385 
386 
387 
389 template<typename Real, typename OtherReal>
390 Real TraceSpSp(const SpMatrix<Real> &A, const SpMatrix<OtherReal> &B);
391 
392 
393 
394 // TraceSpSpLower is the same as Trace(A B) except the lower-diagonal elements
395 // are counted only once not twice as they should be. It is useful in certain
396 // optimizations.
397 template<typename Real>
398 Real TraceSpSpLower(const SpMatrix<Real> &A, const SpMatrix<Real> &B);
399 
400 
403 template<typename Real>
404 Real TraceSpMat(const SpMatrix<Real> &A, const MatrixBase<Real> &B);
405 
408 template<typename Real>
410  const SpMatrix<Real> &B, const MatrixBase<Real> &C,
411  MatrixTransposeType transC);
412 
415 template<typename Real>
417  const SpMatrix<Real> &B, const MatrixBase<Real> &C,
418  MatrixTransposeType transC, const SpMatrix<Real> &D);
419 
424 template<typename Real>
427 Real VecSpVec(const VectorBase<Real> &v1, const SpMatrix<Real> &M,
428  const VectorBase<Real> &v2);
429 
430 
432 
435 
436 
444  BaseFloat K; // maximum condition number
446  std::string name;
450  explicit SolverOptions(const std::string &name):
451  K(1.0e+4), eps(1.0e-40), name(name),
452  optimize_delta(true), diagonal_precondition(false),
453  print_debug_output(true) { }
454  SolverOptions(): K(1.0e+4), eps(1.0e-40), name("[unknown]"),
455  optimize_delta(true), diagonal_precondition(false),
456  print_debug_output(true) { }
457  void Check() const;
458 };
459 
460 
467 
468 template<typename Real>
470  const VectorBase<Real> &g,
471  const SolverOptions &opts,
472  VectorBase<Real> *x);
473 
474 
475 
486 template<typename Real>
488  const MatrixBase<Real> &Y,
489  const SpMatrix<Real> &P,
490  const SolverOptions &opts,
491  MatrixBase<Real> *M);
492 
497 template<typename Real>
499  const SpMatrix<Real> &P1,
500  const SpMatrix<Real> &P2,
501  const SpMatrix<Real> &Q1,
502  const SpMatrix<Real> &Q2,
503  const SolverOptions &opts,
504  MatrixBase<Real> *M);
505 
506 
508 
509 } // namespace kaldi
510 
511 
512 // Including the implementation (now actually just includes some
513 // template specializations).
514 #include "matrix/sp-matrix-inl.h"
515 
516 
517 #endif // KALDI_MATRIX_SP_MATRIX_H_
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
bool IsDiagonal(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:465
Packed symetric matrix class.
Definition: matrix-common.h:62
SpMatrix(const SpMatrix< OtherReal > &orig)
Definition: sp-matrix.h:61
MatrixResizeType
Definition: matrix-common.h:37
This class describes the options for maximizing various quadratic objective functions.
Definition: sp-matrix.h:443
double SolveQuadraticProblem(const SpMatrix< double > &H, const VectorBase< double > &g, const SolverOptions &opts, VectorBase< double > *x)
Definition: sp-matrix.cc:635
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 Qr(MatrixBase< Real > *Q)
The symmetric QR algorithm.
Definition: qr.cc:409
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
Real Cond() const
Returns condition number by computing Svd.
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
Real Trace() const
Definition: sp-matrix.cc:171
SolverOptions(const std::string &name)
Definition: sp-matrix.h:450
SpMatrix(const SpMatrix< Real > &orig)
Definition: sp-matrix.h:57
void InvertDouble(Real *logdet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
Definition: sp-matrix.cc:312
Real FrobeniusNorm() const
sqrt of sum of square elements.
Definition: sp-matrix.cc:513
void AddPacked(const Real alpha, const PackedMatrix< Real > &M)
Real TraceSpMat(const SpMatrix< Real > &A, const MatrixBase< Real > &B)
Returns tr(A B).
Definition: sp-matrix.cc:389
SpMatrix(const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
Definition: sp-matrix.h:71
void swap(basic_filebuf< CharT, Traits > &x, basic_filebuf< CharT, Traits > &y)
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
bool IsTridiagonal(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:491
MatrixIndexT NumRows() const
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
void Tridiagonalize(MatrixBase< Real > *Q)
Tridiagonalize the matrix with an orthogonal transformation.
Definition: qr.cc:165
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 PrintEigs(const char *name)
Definition: sp-matrix.h:203
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 AddVec2(const Real alpha, const VectorBase< OtherReal > &v)
rank-one update, this <– this + alpha v v&#39;
Definition: sp-matrix.cc:946
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 operator()(MatrixIndexT r, MatrixIndexT c) const
Definition: sp-matrix.h:102
Real LogPosDefDet() const
Computes log determinant but only for +ve-def matrices (it uses Cholesky).
Definition: sp-matrix.cc:36
Packed matrix: base class for triangular and symmetric matrices.
Definition: matrix-common.h:64
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
MatrixIndexT LimitCond(Real maxCond=1.0e+5, bool invert=false)
Definition: sp-matrix.cc:606
void AddSp(const Real alpha, const SpMatrix< Real > &Ma)
Definition: sp-matrix.h:211
SpMatrix(MatrixIndexT r, MatrixResizeType resize_type=kSetZero)
Definition: sp-matrix.h:54
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
std::string name
Definition: sp-matrix.h:446
void CopyFromSp(const SpMatrix< OtherReal > &other)
Definition: sp-matrix.h:90
void TopEigs(VectorBase< Real > *s, MatrixBase< Real > *P, MatrixIndexT lanczos_dim=0) const
This function gives you, approximately, the largest eigenvalues of the symmetric matrix and the corre...
Definition: qr.cc:454
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
Real TraceSpSpLower(const SpMatrix< Real > &A, const SpMatrix< Real > &B)
Definition: sp-matrix.cc:1171
Packed symetric matrix class.
Definition: matrix-common.h:63
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 CopyFromMat(const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
Definition: sp-matrix.cc:112
PackedMatrix< Real > & operator=(const PackedMatrix< Real > &other)
Definition: packed-matrix.h:68
void CopyFromPacked(const PackedMatrix< OtherReal > &orig)
A class representing a vector.
Definition: kaldi-vector.h:406
void EigInternal(VectorBase< Real > *s, MatrixBase< Real > *P, Real tolerance, int recurse) 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
Real Cond() const
Returns maximum ratio of singular values.
Definition: sp-matrix.h:145
MatrixTransposeType
Definition: matrix-common.h:32
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
void Resize(MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
Definition: sp-matrix.h:81
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 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
bool IsZero(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:507
SpMatrix< Real > & operator=(const SpMatrix< Real > &other)
Definition: sp-matrix.h:126
Real MaxAbsEig() const
Returns the maximum of the absolute values of any of the eigenvalues.
Definition: sp-matrix.cc:67
void AddDiagVec(const Real alpha, const VectorBase< OtherReal > &v)
diagonal update, this <– this + diag(v)
Definition: sp-matrix.cc:183
Real LogDet(Real *det_sign=NULL) const
Definition: sp-matrix.cc:579
MatrixIndexT LimitCondDouble(Real maxCond=1.0e+5, bool invert=false)
Definition: sp-matrix.h:333
void Resize(MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
Set packed matrix to a specified size (can be zero).