diff options
Diffstat (limited to 'kaldi_io/src/kaldi/matrix/sp-matrix.h')
| -rw-r--r-- | kaldi_io/src/kaldi/matrix/sp-matrix.h | 524 |
1 files changed, 0 insertions, 524 deletions
diff --git a/kaldi_io/src/kaldi/matrix/sp-matrix.h b/kaldi_io/src/kaldi/matrix/sp-matrix.h deleted file mode 100644 index 209d24a..0000000 --- a/kaldi_io/src/kaldi/matrix/sp-matrix.h +++ /dev/null @@ -1,524 +0,0 @@ -// matrix/sp-matrix.h - -// Copyright 2009-2011 Ondrej Glembek; Microsoft Corporation; Lukas Burget; -// Saarland University; Ariya Rastrow; Yanmin Qian; -// Jan Silovsky - -// See ../../COPYING for clarification regarding multiple authors -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at - -// http://www.apache.org/licenses/LICENSE-2.0 - -// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY -// KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED -// WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE, -// MERCHANTABLITY OR NON-INFRINGEMENT. -// See the Apache 2 License for the specific language governing permissions and -// limitations under the License. -#ifndef KALDI_MATRIX_SP_MATRIX_H_ -#define KALDI_MATRIX_SP_MATRIX_H_ - -#include <algorithm> -#include <vector> - -#include "matrix/packed-matrix.h" - -namespace kaldi { - - -/// \addtogroup matrix_group -/// @{ -template<typename Real> class SpMatrix; - - -/** - * @brief Packed symetric matrix class -*/ -template<typename Real> -class SpMatrix : public PackedMatrix<Real> { - friend class CuSpMatrix<Real>; - public: - // so it can use our assignment operator. - friend class std::vector<Matrix<Real> >; - - SpMatrix(): PackedMatrix<Real>() {} - - /// Copy constructor from CUDA version of SpMatrix - /// This is defined in ../cudamatrix/cu-sp-matrix.h - - explicit SpMatrix(const CuSpMatrix<Real> &cu); - - explicit SpMatrix(MatrixIndexT r, MatrixResizeType resize_type = kSetZero) - : PackedMatrix<Real>(r, resize_type) {} - - SpMatrix(const SpMatrix<Real> &orig) - : PackedMatrix<Real>(orig) {} - - template<typename OtherReal> - explicit SpMatrix(const SpMatrix<OtherReal> &orig) - : PackedMatrix<Real>(orig) {} - -#ifdef KALDI_PARANOID - explicit SpMatrix(const MatrixBase<Real> & orig, - SpCopyType copy_type = kTakeMeanAndCheck) - : PackedMatrix<Real>(orig.NumRows(), kUndefined) { - CopyFromMat(orig, copy_type); - } -#else - explicit SpMatrix(const MatrixBase<Real> & orig, - SpCopyType copy_type = kTakeMean) - : PackedMatrix<Real>(orig.NumRows(), kUndefined) { - CopyFromMat(orig, copy_type); - } -#endif - - /// Shallow swap. - void Swap(SpMatrix *other); - - inline void Resize(MatrixIndexT nRows, MatrixResizeType resize_type = kSetZero) { - PackedMatrix<Real>::Resize(nRows, resize_type); - } - - void CopyFromSp(const SpMatrix<Real> &other) { - PackedMatrix<Real>::CopyFromPacked(other); - } - - template<typename OtherReal> - void CopyFromSp(const SpMatrix<OtherReal> &other) { - PackedMatrix<Real>::CopyFromPacked(other); - } - -#ifdef KALDI_PARANOID - void CopyFromMat(const MatrixBase<Real> &orig, - SpCopyType copy_type = kTakeMeanAndCheck); -#else // different default arg if non-paranoid mode. - void CopyFromMat(const MatrixBase<Real> &orig, - SpCopyType copy_type = kTakeMean); -#endif - - inline Real operator() (MatrixIndexT r, MatrixIndexT c) const { - // if column is less than row, then swap these as matrix is stored - // as upper-triangular... only allowed for const matrix object. - if (static_cast<UnsignedMatrixIndexT>(c) > - static_cast<UnsignedMatrixIndexT>(r)) - std::swap(c, r); - // c<=r now so don't have to check c. - KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) < - static_cast<UnsignedMatrixIndexT>(this->num_rows_)); - return *(this->data_ + (r*(r+1)) / 2 + c); - // Duplicating code from PackedMatrix.h - } - - inline Real &operator() (MatrixIndexT r, MatrixIndexT c) { - if (static_cast<UnsignedMatrixIndexT>(c) > - static_cast<UnsignedMatrixIndexT>(r)) - std::swap(c, r); - // c<=r now so don't have to check c. - KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) < - static_cast<UnsignedMatrixIndexT>(this->num_rows_)); - return *(this->data_ + (r * (r + 1)) / 2 + c); - // Duplicating code from PackedMatrix.h - } - - using PackedMatrix<Real>::operator =; - using PackedMatrix<Real>::Scale; - - /// matrix inverse. - /// if inverse_needed = false, will fill matrix with garbage. - /// (only useful if logdet wanted). - void Invert(Real *logdet = NULL, Real *det_sign= NULL, - bool inverse_needed = true); - - // Below routine does inversion in double precision, - // even for single-precision object. - void InvertDouble(Real *logdet = NULL, Real *det_sign = NULL, - bool inverse_needed = true); - - /// Returns maximum ratio of singular values. - inline Real Cond() const { - Matrix<Real> tmp(*this); - return tmp.Cond(); - } - - /// Takes matrix to a fraction power via Svd. - /// Will throw exception if matrix is not positive semidefinite - /// (to within a tolerance) - void ApplyPow(Real exponent); - - /// This is the version of SVD that we implement for symmetric positive - /// definite matrices. This exists for historical reasons; right now its - /// internal implementation is the same as Eig(). It computes the eigenvalue - /// decomposition (*this) = P * diag(s) * P^T with P orthogonal. Will throw - /// exception if input is not positive semidefinite to within a tolerance. - void SymPosSemiDefEig(VectorBase<Real> *s, MatrixBase<Real> *P, - Real tolerance = 0.001) const; - - /// Solves the symmetric eigenvalue problem: at end we should have (*this) = P - /// * diag(s) * P^T. We solve the problem using the symmetric QR method. - /// P may be NULL. - /// Implemented in qr.cc. - /// If you need the eigenvalues sorted, the function SortSvd declared in - /// kaldi-matrix is suitable. - void Eig(VectorBase<Real> *s, MatrixBase<Real> *P = NULL) const; - - /// This function gives you, approximately, the largest eigenvalues of the - /// symmetric matrix and the corresponding eigenvectors. (largest meaning, - /// further from zero). It does this by doing a SVD within the Krylov - /// subspace generated by this matrix and a random vector. This is - /// a form of the Lanczos method with complete reorthogonalization, followed - /// by SVD within a smaller dimension ("lanczos_dim"). - /// - /// If *this is m by m, s should be of dimension n and P should be of - /// dimension m by n, with n <= m. The *columns* of P are the approximate - /// eigenvectors; P * diag(s) * P^T would be a low-rank reconstruction of - /// *this. The columns of P will be orthogonal, and the elements of s will be - /// the eigenvalues of *this projected into that subspace, but beyond that - /// there are no exact guarantees. (This is because the convergence of this - /// method is statistical). Note: it only makes sense to use this - /// method if you are in very high dimension and n is substantially smaller - /// than m: for example, if you want the 100 top eigenvalues of a 10k by 10k - /// matrix. This function calls Rand() to initialize the lanczos - /// iterations and also for restarting. - /// If lanczos_dim is zero, it will default to the greater of: - /// s->Dim() + 50 or s->Dim() + s->Dim()/2, but not more than this->Dim(). - /// If lanczos_dim == this->Dim(), you might as well just call the function - /// Eig() since the result will be the same, and Eig() would be faster; the - /// whole point of this function is to reduce the dimension of the SVD - /// computation. - void TopEigs(VectorBase<Real> *s, MatrixBase<Real> *P, - MatrixIndexT lanczos_dim = 0) const; - - - - /// Takes log of the matrix (does eigenvalue decomposition then takes - /// log of eigenvalues and reconstructs). Will throw of not +ve definite. - void Log(); - - - // Takes exponential of the matrix (equivalent to doing eigenvalue - // decomposition then taking exp of eigenvalues and reconstructing). - void Exp(); - - /// Returns the maximum of the absolute values of any of the - /// eigenvalues. - Real MaxAbsEig() const; - - void PrintEigs(const char *name) { - Vector<Real> s((*this).NumRows()); - Matrix<Real> P((*this).NumRows(), (*this).NumCols()); - SymPosSemiDefEig(&s, &P); - KALDI_LOG << "PrintEigs: " << name << ": " << s; - } - - bool IsPosDef() const; // returns true if Cholesky succeeds. - void AddSp(const Real alpha, const SpMatrix<Real> &Ma) { - this->AddPacked(alpha, Ma); - } - - /// Computes log determinant but only for +ve-def matrices - /// (it uses Cholesky). - /// If matrix is not +ve-def, it will throw an exception - /// was LogPDDeterminant() - Real LogPosDefDet() const; - - Real LogDet(Real *det_sign = NULL) const; - - /// rank-one update, this <-- this + alpha v v' - template<typename OtherReal> - void AddVec2(const Real alpha, const VectorBase<OtherReal> &v); - - /// rank-two update, this <-- this + alpha (v w' + w v'). - void AddVecVec(const Real alpha, const VectorBase<Real> &v, - const VectorBase<Real> &w); - - /// Does *this = beta * *thi + alpha * diag(v) * S * diag(v) - void AddVec2Sp(const Real alpha, const VectorBase<Real> &v, - const SpMatrix<Real> &S, const Real beta); - - /// diagonal update, this <-- this + diag(v) - template<typename OtherReal> - void AddDiagVec(const Real alpha, const VectorBase<OtherReal> &v); - - /// rank-N update: - /// if (transM == kNoTrans) - /// (*this) = beta*(*this) + alpha * M * M^T, - /// or (if transM == kTrans) - /// (*this) = beta*(*this) + alpha * M^T * M - /// Note: beta used to default to 0.0. - void AddMat2(const Real alpha, const MatrixBase<Real> &M, - MatrixTransposeType transM, const Real beta); - - /// Extension of rank-N update: - /// this <-- beta*this + alpha * M * A * M^T. - /// (*this) and A are allowed to be the same. - /// If transM == kTrans, then we do it as M^T * A * M. - void AddMat2Sp(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. - /// This was required for making the raw-fMLLR code efficient. - void AddSmat2Sp(const Real alpha, const MatrixBase<Real> &M, - MatrixTransposeType transM, const SpMatrix<Real> &A, - const Real beta = 0.0); - - /// The following function does: - /// this <-- beta*this + alpha * T * A * T^T. - /// (*this) and A are allowed to be the same. - /// If transM == kTrans, then we do it as alpha * T^T * A * T. - /// Currently it just calls AddMat2Sp, but if needed we - /// can implement it more efficiently. - 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 * T^T. - /// (*this) and A are allowed to be the same. - /// If transM == kTrans, then we do it as alpha * T^T * T - /// Currently it just calls AddMat2, but if needed we - /// can implement it more efficiently. - void AddTp2(const Real alpha, const TpMatrix<Real> &T, - MatrixTransposeType transM, const Real beta = 0.0); - - /// Extension of rank-N update: - /// this <-- beta*this + alpha * M * diag(v) * M^T. - /// if transM == kTrans, then - /// this <-- beta*this + alpha * M^T * diag(v) * M. - void AddMat2Vec(const Real alpha, const MatrixBase<Real> &M, - MatrixTransposeType transM, const VectorBase<Real> &v, - const Real beta = 0.0); - - - /// Floors this symmetric matrix to the matrix - /// alpha * Floor, where the matrix Floor is positive - /// definite. - /// It is floored in the sense that after flooring, - /// x^T (*this) x >= x^T (alpha*Floor) x. - /// This is accomplished using an Svd. It will crash - /// if Floor is not positive definite. Returns the number of - /// elements that were floored. - int ApplyFloor(const SpMatrix<Real> &Floor, Real alpha = 1.0, - bool verbose = false); - - /// Floor: Given a positive semidefinite matrix, floors the eigenvalues - /// to the specified quantity. A previous version of this function had - /// a tolerance which is now no longer needed since we have code to - /// do the symmetric eigenvalue decomposition and no longer use the SVD - /// code for that purose. - int ApplyFloor(Real floor); - - bool IsDiagonal(Real cutoff = 1.0e-05) const; - bool IsUnit(Real cutoff = 1.0e-05) const; - bool IsZero(Real cutoff = 1.0e-05) const; - bool IsTridiagonal(Real cutoff = 1.0e-05) const; - - /// sqrt of sum of square elements. - Real FrobeniusNorm() const; - - /// Returns true if ((*this)-other).FrobeniusNorm() <= - /// tol*(*this).FrobeniusNorma() - bool ApproxEqual(const SpMatrix<Real> &other, float tol = 0.01) const; - - // LimitCond: - // Limits the condition of symmetric positive semidefinite matrix to - // a specified value - // by flooring all eigenvalues to a positive number which is some multiple - // of the largest one (or zero if there are no positive eigenvalues). - // Takes the condition number we are willing to accept, and floors - // eigenvalues to the largest eigenvalue divided by this. - // Returns #eigs floored or already equal to the floor. - // Throws exception if input is not positive definite. - // returns #floored. - MatrixIndexT LimitCond(Real maxCond = 1.0e+5, bool invert = false); - - // as LimitCond but all done in double precision. // returns #floored. - MatrixIndexT LimitCondDouble(Real maxCond = 1.0e+5, bool invert = false) { - SpMatrix<double> dmat(*this); - MatrixIndexT ans = dmat.LimitCond(maxCond, invert); - (*this).CopyFromSp(dmat); - return ans; - } - Real Trace() const; - - /// Tridiagonalize the matrix with an orthogonal transformation. If - /// *this starts as S, produce T (and Q, if non-NULL) such that - /// T = Q A Q^T, i.e. S = Q^T T Q. Caution: this is the other way - /// round from most authors (it's more efficient in row-major indexing). - void Tridiagonalize(MatrixBase<Real> *Q); - - /// The symmetric QR algorithm. This will mostly be useful in internal code. - /// Typically, you will call this after Tridiagonalize(), on the same object. - /// When called, *this (call it A at this point) must be tridiagonal; at exit, - /// *this will be a diagonal matrix D that is similar to A via orthogonal - /// transformations. This algorithm right-multiplies Q by orthogonal - /// transformations. It turns *this from a tridiagonal into a diagonal matrix - /// while maintaining that (Q *this Q^T) has the same value at entry and exit. - /// At entry Q should probably be either NULL or orthogonal, but we don't check - /// this. - void Qr(MatrixBase<Real> *Q); - - private: - void EigInternal(VectorBase<Real> *s, MatrixBase<Real> *P, - Real tolerance, int recurse) const; -}; - -/// @} end of "addtogroup matrix_group" - -/// \addtogroup matrix_funcs_scalar -/// @{ - - -/// Returns tr(A B). -float TraceSpSp(const SpMatrix<float> &A, const SpMatrix<float> &B); -double TraceSpSp(const SpMatrix<double> &A, const SpMatrix<double> &B); - - -template<typename Real> -inline bool ApproxEqual(const SpMatrix<Real> &A, - const SpMatrix<Real> &B, Real tol = 0.01) { - return A.ApproxEqual(B, tol); -} - -template<typename Real> -inline void AssertEqual(const SpMatrix<Real> &A, - const SpMatrix<Real> &B, Real tol = 0.01) { - KALDI_ASSERT(ApproxEqual(A, B, tol)); -} - - - -/// Returns tr(A B). -template<typename Real, typename OtherReal> -Real TraceSpSp(const SpMatrix<Real> &A, const SpMatrix<OtherReal> &B); - - - -// TraceSpSpLower is the same as Trace(A B) except the lower-diagonal elements -// are counted only once not twice as they should be. It is useful in certain -// optimizations. -template<typename Real> -Real TraceSpSpLower(const SpMatrix<Real> &A, const SpMatrix<Real> &B); - - -/// Returns tr(A B). -/// No option to transpose B because would make no difference. -template<typename Real> -Real TraceSpMat(const SpMatrix<Real> &A, const MatrixBase<Real> &B); - -/// Returns tr(A B C) -/// (A and C may be transposed as specified by transA and transC). -template<typename Real> -Real TraceMatSpMat(const MatrixBase<Real> &A, MatrixTransposeType transA, - const SpMatrix<Real> &B, const MatrixBase<Real> &C, - MatrixTransposeType transC); - -/// Returns tr (A B C D) -/// (A and C may be transposed as specified by transA and transB). -template<typename Real> -Real TraceMatSpMatSp(const MatrixBase<Real> &A, MatrixTransposeType transA, - const SpMatrix<Real> &B, const MatrixBase<Real> &C, - MatrixTransposeType transC, const SpMatrix<Real> &D); - -/** Computes v1^T * M * v2. Not as efficient as it could be where v1 == v2 - * (but no suitable blas routines available). - */ - -/// Returns \f$ v_1^T M v_2 \f$ -/// Not as efficient as it could be where v1 == v2. -template<typename Real> -Real VecSpVec(const VectorBase<Real> &v1, const SpMatrix<Real> &M, - const VectorBase<Real> &v2); - - -/// @} \addtogroup matrix_funcs_scalar - -/// \addtogroup matrix_funcs_misc -/// @{ - - -/// This class describes the options for maximizing various quadratic objective -/// functions. It's mostly as described in the SGMM paper "the subspace -/// Gaussian mixture model -- a structured model for speech recognition", but -/// the diagonal_precondition option is newly added, to handle problems where -/// different dimensions have very different scaling (we recommend to use the -/// option but it's set false for back compatibility). -struct SolverOptions { - BaseFloat K; // maximum condition number - BaseFloat eps; - std::string name; - bool optimize_delta; - bool diagonal_precondition; - bool print_debug_output; - explicit SolverOptions(const std::string &name): - K(1.0e+4), eps(1.0e-40), name(name), - optimize_delta(true), diagonal_precondition(false), - print_debug_output(true) { } - SolverOptions(): K(1.0e+4), eps(1.0e-40), name("[unknown]"), - optimize_delta(true), diagonal_precondition(false), - print_debug_output(true) { } - void Check() const; -}; - - -/// Maximizes the auxiliary function -/// \f[ Q(x) = x.g - 0.5 x^T H x \f] -/// using a numerically stable method. Like a numerically stable version of -/// \f$ x := Q^{-1} g. \f$ -/// Assumes H positive semidefinite. -/// Returns the objective-function change. - -template<typename Real> -Real SolveQuadraticProblem(const SpMatrix<Real> &H, - const VectorBase<Real> &g, - const SolverOptions &opts, - VectorBase<Real> *x); - - - -/// Maximizes the auxiliary function : -/// \f[ Q(x) = tr(M^T P Y) - 0.5 tr(P M Q M^T) \f] -/// Like a numerically stable version of \f$ M := Y Q^{-1} \f$. -/// Assumes Q and P positive semidefinite, and matrix dimensions match -/// enough to make expressions meaningful. -/// This is mostly as described in the SGMM paper "the subspace Gaussian mixture -/// model -- a structured model for speech recognition", but the -/// diagonal_precondition option is newly added, to handle problems -/// where different dimensions have very different scaling (we recommend to use -/// the option but it's set false for back compatibility). -template<typename Real> -Real SolveQuadraticMatrixProblem(const SpMatrix<Real> &Q, - const MatrixBase<Real> &Y, - const SpMatrix<Real> &P, - const SolverOptions &opts, - MatrixBase<Real> *M); - -/// Maximizes the auxiliary function : -/// \f[ Q(M) = tr(M^T G) -0.5 tr(P_1 M Q_1 M^T) -0.5 tr(P_2 M Q_2 M^T). \f] -/// Encountered in matrix update with a prior. We also apply a limit on the -/// condition but it should be less frequently necessary, and can be set larger. -template<typename Real> -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); - - -/// @} End of "addtogroup matrix_funcs_misc" - -} // namespace kaldi - - -// Including the implementation (now actually just includes some -// template specializations). -#include "matrix/sp-matrix-inl.h" - - -#endif // KALDI_MATRIX_SP_MATRIX_H_ - |