All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
SpMatrix< Real > Class Template Reference

Packed symetric matrix class. More...

#include <matrix-common.h>

Inheritance diagram for SpMatrix< Real >:
Collaboration diagram for SpMatrix< Real >:

Public Member Functions

 SpMatrix ()
 
 SpMatrix (const CuSpMatrix< Real > &cu)
 Copy constructor from CUDA version of SpMatrix This is defined in ../cudamatrix/cu-sp-matrix.h. More...
 
 SpMatrix (MatrixIndexT r, MatrixResizeType resize_type=kSetZero)
 
 SpMatrix (const SpMatrix< Real > &orig)
 
template<typename OtherReal >
 SpMatrix (const SpMatrix< OtherReal > &orig)
 
 SpMatrix (const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
 
void Swap (SpMatrix *other)
 Shallow swap. More...
 
void Resize (MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
 
void CopyFromSp (const SpMatrix< Real > &other)
 
template<typename OtherReal >
void CopyFromSp (const SpMatrix< OtherReal > &other)
 
void CopyFromMat (const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
 
Real operator() (MatrixIndexT r, MatrixIndexT c) const
 
Real & operator() (MatrixIndexT r, MatrixIndexT c)
 
SpMatrix< Real > & operator= (const SpMatrix< Real > &other)
 
void Invert (Real *logdet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
 matrix inverse. More...
 
void InvertDouble (Real *logdet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
 
Real Cond () const
 Returns maximum ratio of singular values. More...
 
void ApplyPow (Real exponent)
 Takes matrix to a fraction power via Svd. More...
 
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. More...
 
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. More...
 
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 corresponding eigenvectors. More...
 
void Log ()
 Takes log of the matrix (does eigenvalue decomposition then takes log of eigenvalues and reconstructs). More...
 
void Exp ()
 
Real MaxAbsEig () const
 Returns the maximum of the absolute values of any of the eigenvalues. More...
 
void PrintEigs (const char *name)
 
bool IsPosDef () const
 
void AddSp (const Real alpha, const SpMatrix< Real > &Ma)
 
Real LogPosDefDet () const
 Computes log determinant but only for +ve-def matrices (it uses Cholesky). More...
 
Real LogDet (Real *det_sign=NULL) const
 
template<typename OtherReal >
void AddVec2 (const Real alpha, const VectorBase< OtherReal > &v)
 rank-one update, this <– this + alpha v v' More...
 
void AddVecVec (const Real alpha, const VectorBase< Real > &v, const VectorBase< Real > &w)
 rank-two update, this <– this + alpha (v w' + w v'). More...
 
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) More...
 
template<typename OtherReal >
void AddDiagVec (const Real alpha, const VectorBase< OtherReal > &v)
 diagonal update, this <– this + diag(v) More...
 
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 == kTrans) (*this) = beta*(*this) + alpha * M^T * M Note: beta used to default to 0.0. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
int ApplyFloor (Real floor)
 Floor: Given a positive semidefinite matrix, floors the eigenvalues to the specified quantity. More...
 
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
 
Real FrobeniusNorm () const
 sqrt of sum of square elements. More...
 
bool ApproxEqual (const SpMatrix< Real > &other, float tol=0.01) const
 Returns true if ((*this)-other).FrobeniusNorm() <= tol*(*this).FrobeniusNorma() More...
 
MatrixIndexT LimitCond (Real maxCond=1.0e+5, bool invert=false)
 
MatrixIndexT LimitCondDouble (Real maxCond=1.0e+5, bool invert=false)
 
Real Trace () const
 
void Tridiagonalize (MatrixBase< Real > *Q)
 Tridiagonalize the matrix with an orthogonal transformation. More...
 
void Qr (MatrixBase< Real > *Q)
 The symmetric QR algorithm. More...
 
template<>
void AddVec2 (const double alpha, const VectorBase< double > &v)
 
template<>
void AddVec2 (const float alpha, const VectorBase< float > &v)
 
template<>
void AddVec2 (const double alpha, const VectorBase< double > &v)
 
- Public Member Functions inherited from PackedMatrix< Real >
 PackedMatrix ()
 
 PackedMatrix (MatrixIndexT r, MatrixResizeType resize_type=kSetZero)
 
 PackedMatrix (const PackedMatrix< Real > &orig)
 
template<typename OtherReal >
 PackedMatrix (const PackedMatrix< OtherReal > &orig)
 
void SetZero ()
 
void SetUnit ()
 < Set to zero More...
 
void SetRandn ()
 < Set to unit matrix. More...
 
Real Trace () const
 < Set to random values of a normal distribution More...
 
PackedMatrix< Real > & operator= (const PackedMatrix< Real > &other)
 
 ~PackedMatrix ()
 
void Resize (MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
 Set packed matrix to a specified size (can be zero). More...
 
void AddToDiag (const Real r)
 
void ScaleDiag (const Real alpha)
 
void SetDiag (const Real alpha)
 
template<typename OtherReal >
void CopyFromPacked (const PackedMatrix< OtherReal > &orig)
 
template<typename OtherReal >
void CopyFromVec (const SubVector< OtherReal > &orig)
 CopyFromVec just interprets the vector as having the same layout as the packed matrix. More...
 
Real * Data ()
 
const Real * Data () const
 
MatrixIndexT NumRows () const
 
MatrixIndexT NumCols () const
 
size_t SizeInBytes () const
 
Real operator() (MatrixIndexT r, MatrixIndexT c) const
 
Real & operator() (MatrixIndexT r, MatrixIndexT c)
 
Real Max () const
 
Real Min () const
 
void Scale (Real c)
 
void Read (std::istream &in, bool binary, bool add=false)
 
void Write (std::ostream &out, bool binary) const
 
void Destroy ()
 
void Swap (PackedMatrix< Real > *other)
 Swaps the contents of *this and *other. Shallow swap. More...
 
void Swap (Matrix< Real > *other)
 

Private Member Functions

void EigInternal (VectorBase< Real > *s, MatrixBase< Real > *P, Real tolerance, int recurse) const
 

Friends

class CuSpMatrix< Real >
 
class std::vector< Matrix< Real > >
 

Additional Inherited Members

- Protected Member Functions inherited from PackedMatrix< Real >
void AddPacked (const Real alpha, const PackedMatrix< Real > &M)
 
- Protected Attributes inherited from PackedMatrix< Real >
Real * data_
 
MatrixIndexT num_rows_
 

Detailed Description

template<typename Real>
class kaldi::SpMatrix< Real >

Packed symetric matrix class.

Definition at line 60 of file matrix-common.h.

Constructor & Destructor Documentation

SpMatrix ( )
inline

Definition at line 47 of file sp-matrix.h.

47 : PackedMatrix<Real>() {}
SpMatrix ( const CuSpMatrix< Real > &  cu)
explicit

Copy constructor from CUDA version of SpMatrix This is defined in ../cudamatrix/cu-sp-matrix.h.

Definition at line 154 of file cu-sp-matrix.h.

154  {
155  Resize(cu.NumRows());
156  cu.CopyToSp(this);
157 }
void Resize(MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
Definition: sp-matrix.h:81
SpMatrix ( MatrixIndexT  r,
MatrixResizeType  resize_type = kSetZero 
)
inlineexplicit

Definition at line 54 of file sp-matrix.h.

55  : PackedMatrix<Real>(r, resize_type) {}
SpMatrix ( const SpMatrix< Real > &  orig)
inline

Definition at line 57 of file sp-matrix.h.

58  : PackedMatrix<Real>(orig) {}
SpMatrix ( const SpMatrix< OtherReal > &  orig)
inlineexplicit

Definition at line 61 of file sp-matrix.h.

62  : PackedMatrix<Real>(orig) {}
SpMatrix ( const MatrixBase< Real > &  orig,
SpCopyType  copy_type = kTakeMean 
)
inlineexplicit

Definition at line 71 of file sp-matrix.h.

73  : PackedMatrix<Real>(orig.NumRows(), kUndefined) {
74  CopyFromMat(orig, copy_type);
75  }
void CopyFromMat(const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
Definition: sp-matrix.cc:138

Member Function Documentation

template void AddDiagVec ( const Real  alpha,
const VectorBase< OtherReal > &  v 
)

diagonal update, this <– this + diag(v)

Definition at line 209 of file sp-matrix.cc.

References VectorBase< Real >::Data(), data_, VectorBase< Real >::Dim(), rnnlm::i, and KALDI_ASSERT.

Referenced by Fmpe::ComputeC(), main(), and kaldi::UnitTestSpAddDiagVec().

209  {
210  int32 num_rows = this->num_rows_;
211  KALDI_ASSERT(num_rows == v.Dim() && num_rows > 0);
212  const OtherReal *src = v.Data();
213  Real *dst = this->data_;
214  if (alpha == 1.0)
215  for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
216  *dst += *src;
217  else
218  for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
219  *dst += alpha * *src;
220 }
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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 == kTrans) (*this) = beta*(*this) + alpha * M^T * M Note: beta used to default to 0.0.

Definition at line 1134 of file sp-matrix.cc.

References kaldi::cblas_Xsyrk(), MatrixBase< Real >::Data(), KALDI_ASSERT, kaldi::kNoTrans, kaldi::kTakeLower, kaldi::kTrans, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), and MatrixBase< Real >::Stride().

Referenced by CovarianceStats::AccStats(), PldaStats::AddSamples(), SpMatrix< Real >::ApplyFloor(), kaldi::CholeskyUnitTestTr(), kaldi::ComputePca(), main(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), rand_posdef_spmatrix(), kaldi::RandFullCova(), kaldi::unittest::RandPosdefSpMatrix(), RandPosdefSpMatrix(), kaldi::RandPosdefSpMatrix(), MatrixBase< Real >::SymPosSemiDefEig(), kaldi::UnitTestAddMat2(), kaldi::UnitTestCuSpMatrixAddMat2(), kaldi::UnitTestCuSpMatrixInvert(), kaldi::UnitTestCuTpMatrixCholesky(), kaldi::UnitTestFloorChol(), kaldi::UnitTestFloorUnit(), kaldi::UnitTestGpsr(), kaldi::UnitTestLimitCondInvert(), kaldi::UnitTestMat2Vec(), kaldi::UnitTestPca(), kaldi::UnitTestPca2(), kaldi::UnitTestPldaEstimation(), kaldi::UnitTestRankNUpdate(), kaldi::UnitTestSolve(), kaldi::UnitTestSpLogExp(), kaldi::UnitTestTopEigs(), and kaldi::UnitTestTp2().

1135  {
1136  KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows())
1137  || (transM == kTrans && this->NumRows() == M.NumCols()));
1138 
1139  // Cblas has no function *sprk (i.e. symmetric packed rank-k update), so we
1140  // use as temporary storage a regular matrix of which we only access its lower
1141  // triangle
1142 
1143  MatrixIndexT this_dim = this->NumRows(),
1144  m_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows());
1145 
1146  if (this_dim == 0) return;
1147  if (alpha == 0.0) {
1148  if (beta != 1.0) this->Scale(beta);
1149  return;
1150  }
1151 
1152  Matrix<Real> temp_mat(*this); // wastefully copies upper triangle too, but this
1153  // doesn't dominate O(N) time.
1154 
1155  // This function call is hard-coded to update the lower triangle.
1156  cblas_Xsyrk(transM, this_dim, m_other_dim, alpha, M.Data(),
1157  M.Stride(), beta, temp_mat.Data(), temp_mat.Stride());
1158 
1159  this->CopyFromMat(temp_mat, kTakeLower);
1160 }
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)
void Scale(Real c)
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
void CopyFromMat(const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
Definition: sp-matrix.cc:138
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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.

(*this) and A are allowed to be the same. If transM == kTrans, then we do it as M^T * A * M.

Definition at line 1006 of file sp-matrix.cc.

References kaldi::cblas_Xgemv(), kaldi::cblas_Xspmv(), SpMatrix< Real >::CopyFromSp(), VectorBase< Real >::Data(), MatrixBase< Real >::Data(), PackedMatrix< Real >::Data(), KALDI_ASSERT, kaldi::kNoTrans, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), PackedMatrix< Real >::NumRows(), SpMatrix< Real >::Resize(), MatrixBase< Real >::RowData(), PackedMatrix< Real >::SizeInBytes(), and MatrixBase< Real >::Stride().

Referenced by kaldi::ApplyFeatureTransformToStats(), SpMatrix< Real >::ApplyFloor(), IvectorExtractor::ComputeDerivedVars(), kaldi::ComputeFeatureNormalizingTransform(), AmSgmm2::ComputeFmllrPreXform(), AmSgmm2::ComputeH(), Sgmm2Project::ComputeLdaTransform(), kaldi::ComputeLdaTransform(), MleAmSgmm2Updater::ComputeMPrior(), Sgmm2Project::ComputeProjection(), LdaEstimate::Estimate(), BasisFmllrEstimate::EstimateFmllrBasis(), FeatureTransformEstimate::EstimateInternal(), FeatureTransformEstimateMulti::EstimateTransformPart(), PldaEstimator::GetOutput(), Sgmm2Project::ProjectVariance(), MleAmSgmm2Updater::RenormalizeV(), kaldi::SolveDoubleQuadraticMatrixProblem(), kaldi::SolveQuadraticMatrixProblem(), kaldi::UnitTestAddVec2Sp(), kaldi::UnitTestPca2(), kaldi::UnitTestTp2Sp(), kaldi::UnitTestTransposeScatter(), kaldi::UnitTestTridiagonalize(), kaldi::UnitTestTridiagonalizeAndQr(), PldaUnsupervisedAdaptor::UpdatePlda(), IvectorExtractorStats::UpdatePrior(), and IvectorExtractorStats::UpdateVariances().

1008  {
1009  if (transM == kNoTrans) {
1010  KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
1011  } else {
1012  KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
1013  }
1014  Vector<Real> tmp_vec(A.NumRows());
1015  Real *tmp_vec_data = tmp_vec.Data();
1016  SpMatrix<Real> tmp_A;
1017  const Real *p_A_data = A.Data();
1018  Real *p_row_data = this->Data();
1019  MatrixIndexT M_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows()),
1020  M_same_dim = (transM == kNoTrans ? M.NumRows() : M.NumCols()),
1021  M_stride = M.Stride(), dim = this->NumRows();
1022  KALDI_ASSERT(M_same_dim == dim);
1023 
1024  const Real *M_data = M.Data();
1025 
1026  if (this->Data() <= A.Data() + A.SizeInBytes() &&
1027  this->Data() + this->SizeInBytes() >= A.Data()) {
1028  // Matrices A and *this overlap. Make copy of A
1029  tmp_A.Resize(A.NumRows());
1030  tmp_A.CopyFromSp(A);
1031  p_A_data = tmp_A.Data();
1032  }
1033 
1034  if (transM == kNoTrans) {
1035  for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
1036  cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.RowData(r), 1, 0.0, tmp_vec_data, 1);
1037  cblas_Xgemv(transM, r+1, M_other_dim, alpha, M_data, M_stride,
1038  tmp_vec_data, 1, beta, p_row_data, 1);
1039  }
1040  } else {
1041  for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
1042  cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.Data() + r, M.Stride(), 0.0, tmp_vec_data, 1);
1043  cblas_Xgemv(transM, M_other_dim, r+1, alpha, M_data, M_stride,
1044  tmp_vec_data, 1, beta, p_row_data, 1);
1045  }
1046  }
1047 }
size_t SizeInBytes() const
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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)
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 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.

if transM == kTrans, then this <– beta*this + alpha * M^T * diag(v) * M.

Definition at line 1105 of file sp-matrix.cc.

References kaldi::cblas_Xspr(), VectorBase< Real >::Data(), MatrixBase< Real >::Data(), data_, VectorBase< Real >::Dim(), KALDI_ASSERT, kaldi::kNoTrans, kaldi::kTrans, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), and MatrixBase< Real >::Stride().

Referenced by IvectorExtractor::GetAcousticAuxfWeight(), IvectorExtractor::GetIvectorDistWeight(), IvectorExtractor::InvertWithFlooring(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), kaldi::UnitTestEigSp(), and kaldi::UnitTestMat2Vec().

1109  {
1110  this->Scale(beta);
1111  KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows() &&
1112  M.NumCols() == v.Dim()) ||
1113  (transM == kTrans && this->NumRows() == M.NumCols() &&
1114  M.NumRows() == v.Dim()));
1115 
1116  if (transM == kNoTrans) {
1117  const Real *Mdata = M.Data(), *vdata = v.Data();
1118  Real *data = this->data_;
1119  MatrixIndexT dim = this->NumRows(), mcols = M.NumCols(),
1120  mstride = M.Stride();
1121  for (MatrixIndexT col = 0; col < mcols; col++, vdata++, Mdata += 1)
1122  cblas_Xspr(dim, *vdata*alpha, Mdata, mstride, data);
1123  } else {
1124  const Real *Mdata = M.Data(), *vdata = v.Data();
1125  Real *data = this->data_;
1126  MatrixIndexT dim = this->NumRows(), mrows = M.NumRows(),
1127  mstride = M.Stride();
1128  for (MatrixIndexT row = 0; row < mrows; row++, vdata++, Mdata += mstride)
1129  cblas_Xspr(dim, *vdata*alpha, Mdata, 1, data);
1130  }
1131 }
void Scale(Real c)
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
void cblas_Xspr(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, float *Adata)
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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.

This was required for making the raw-fMLLR code efficient.

Definition at line 1050 of file sp-matrix.cc.

References MatrixBase< Real >::AddSmatMat(), MatrixBase< Real >::Data(), data_, rnnlm::i, KALDI_ASSERT, kaldi::kNoTrans, kaldi::kTrans, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), PackedMatrix< Real >::NumRows(), MatrixBase< Real >::Stride(), and kaldi::Xgemv_sparsevec().

Referenced by FmllrRawAccs::ConvertToSimpleStats(), and kaldi::UnitTestAddMat2Sp().

1053  {
1054  KALDI_ASSERT((transM == kNoTrans && M.NumCols() == A.NumRows()) ||
1055  (transM == kTrans && M.NumRows() == A.NumRows()));
1056  if (transM == kNoTrans) {
1057  KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
1058  } else {
1059  KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
1060  }
1061  MatrixIndexT Adim = A.NumRows(), dim = this->num_rows_;
1062 
1063  Matrix<Real> temp_A(A); // represent A as full matrix.
1064  Matrix<Real> temp_MA(dim, Adim);
1065  temp_MA.AddSmatMat(1.0, M, transM, temp_A, kNoTrans, 0.0);
1066 
1067  // Next-- we want to do *this = alpha * temp_MA * M^T + beta * *this.
1068  // To make it sparse vector multiplies, since M is sparse, we'd like
1069  // to do: for each column c, (*this column c) += temp_MA * (M^T's column c.)
1070  // [ignoring the alpha and beta here.]
1071  // It's not convenient to process columns in the symmetric
1072  // packed format because they don't have a constant stride. However,
1073  // we can use the fact that temp_MA * M is symmetric, to just assign
1074  // each row of *this instead of each column.
1075  // So the final iteration is:
1076  // for i = 0... dim-1,
1077  // [the i'th row of *this] = beta * [the i'th row of *this] + alpha *
1078  // temp_MA * [the i'th column of M].
1079  // Of course, we only process the first 0 ... i elements of this row,
1080  // as that's all that are kept in the symmetric packed format.
1081 
1082  Matrix<Real> temp_this(*this);
1083  Real *data = this->data_;
1084  const Real *Mdata = M.Data(), *MAdata = temp_MA.Data();
1085  MatrixIndexT temp_MA_stride = temp_MA.Stride(), Mstride = M.Stride();
1086 
1087  if (transM == kNoTrans) {
1088  // The column of M^T corresponds to the rows of the supplied matrix.
1089  for (MatrixIndexT i = 0; i < dim; i++, data += i) {
1090  MatrixIndexT num_rows = i + 1, num_cols = Adim;
1091  Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
1092  temp_MA_stride, Mdata + (i * Mstride), 1, beta, data, 1);
1093  }
1094  } else {
1095  // The column of M^T corresponds to the columns of the supplied matrix.
1096  for (MatrixIndexT i = 0; i < dim; i++, data += i) {
1097  MatrixIndexT num_rows = i + 1, num_cols = Adim;
1098  Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
1099  temp_MA_stride, Mdata + i, Mstride, beta, data, 1);
1100  }
1101  }
1102 }
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)
int32 MatrixIndexT
Definition: matrix-common.h:96
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void AddSp ( const Real  alpha,
const SpMatrix< Real > &  Ma 
)
inline

Definition at line 221 of file sp-matrix.h.

Referenced by IvectorExtractorStats::Add(), CovarianceStats::AddStats(), SpMatrix< Real >::ApproxEqual(), IvectorExtractorStats::CommitStatsForPrior(), kaldi::ComputeFeatureNormalizingTransform(), AmSgmm2::ComputeHsmFromModel(), Sgmm2Project::ComputeLdaStats(), kaldi::ComputeLdaTransform(), PldaEstimator::ComputeObjfPart2(), EbwAmSgmm2Updater::ComputePhoneVecStats(), LdaEstimate::Estimate(), FeatureTransformEstimate::EstimateInternal(), PldaEstimator::GetStatsFromClassMeans(), PldaEstimator::GetStatsFromIntraClass(), CovarianceStats::GetWithinCovar(), FullGmm::MergedComponentsLogdet(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), LdaEstimate::Read(), MleAmSgmm2Updater::RenormalizeN(), MleAmSgmm2Updater::RenormalizeV(), kaldi::SolveDoubleQuadraticMatrixProblem(), kaldi::UnitTestCuSpMatrixAddSp(), kaldi::UnitTestCuSpMatrixApproxEqual(), kaldi::UnitTestEigSp(), UnitTestEstimateLda(), EbwAmSgmm2Updater::UpdateM(), EbwAmSgmm2Updater::UpdateN(), MleSgmm2SpeakerAccs::UpdateNoU(), EbwAmSgmm2Updater::UpdatePhoneVectorsInternal(), MleAmSgmm2Updater::UpdatePhoneVectorsInternal(), EbwAmSgmm2Updater::UpdateU(), IvectorExtractorStats::UpdateVariances(), EbwAmSgmm2Updater::UpdateVars(), MleAmSgmm2Updater::UpdateVars(), MleSgmm2SpeakerAccs::UpdateWithU(), and FisherComputationClass::~FisherComputationClass().

221  {
222  this->AddPacked(alpha, Ma);
223  }
void AddPacked(const Real alpha, const PackedMatrix< Real > &M)
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.

(*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.

Definition at line 1180 of file sp-matrix.cc.

Referenced by CompressedAffineXformStats::ExtractOneG(), main(), and CuMatrixBase< Real >::SymInvertPosDef().

1181  {
1182  Matrix<Real> Tmat(T);
1183  AddMat2(alpha, Tmat, transM, beta);
1184 }
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:1134
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.

(*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.

Definition at line 1163 of file sp-matrix.cc.

Referenced by Sgmm2Project::ComputeLdaTransform(), kaldi::UnitTestPldaEstimation(), and PldaUnsupervisedAdaptor::UpdatePlda().

1165  {
1166  Matrix<Real> Tmat(T);
1167  AddMat2Sp(alpha, Tmat, transM, A, beta);
1168 }
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:1006
void AddVec2 ( const Real  alpha,
const VectorBase< OtherReal > &  v 
)

rank-one update, this <– this + alpha v v'

Definition at line 970 of file sp-matrix.cc.

References VectorBase< Real >::Data(), data_, VectorBase< Real >::Dim(), rnnlm::i, rnnlm::j, and KALDI_ASSERT.

Referenced by CovarianceStats::AccStats(), IvectorExtractorUtteranceStats::AccStats(), BasisFmllrAccus::AccuGradientScatter(), LdaEstimate::Accumulate(), FmllrSgmm2Accs::AccumulateForFmllrSubspace(), RegtreeMllrDiagGmmAccs::AccumulateForGaussian(), RegtreeFmllrDiagGmmAccs::AccumulateForGaussian(), RegtreeMllrDiagGmmAccs::AccumulateForGmm(), RegtreeFmllrDiagGmmAccs::AccumulateForGmm(), AccumFullGmm::AccumulateFromPosteriors(), FmllrSgmm2Accs::AccumulateFromPosteriors(), PldaStats::AddSamples(), PldaUnsupervisedAdaptor::AddStats(), FmllrRawAccs::CommitSingleFrameStats(), FmllrDiagGmmAccs::CommitSingleFrameStats(), IvectorExtractorStats::CommitStatsForM(), IvectorExtractorStats::CommitStatsForPrior(), IvectorExtractorStats::CommitStatsForWPoint(), Fmpe::ComputeC(), kaldi::ComputeFeatureNormalizingTransform(), AmSgmm2::ComputeFmllrPreXform(), Sgmm2Project::ComputeLdaStats(), EbwAmSgmm2Updater::ComputePhoneVecStats(), IvectorExtractor::GetAcousticAuxfVariance(), LdaEstimate::GetStats(), PldaEstimator::GetStatsFromClassMeans(), FullGmm::LogLikelihoodsPreselect(), main(), FullGmm::MergedComponentsLogdet(), FisherComputationClass::operator()(), IvectorExtractorStats::PriorDiagnostics(), AccumFullGmm::Read(), LdaEstimate::Read(), MleAmSgmm2Updater::RenormalizeV(), kaldi::UnitTestCuSpMatrixAddVec2(), UnitTestEstimateLda(), kaldi::UnitTestSimpleForMat(), MleAmSgmm2Updater::UpdatePhoneVectorsInternal(), PldaUnsupervisedAdaptor::UpdatePlda(), IvectorExtractorStats::UpdatePrior(), MleAmSgmm2Updater::UpdateWGetStats(), MleSgmm2SpeakerAccs::UpdateWithU(), AccumFullGmm::Write(), and LdaEstimate::Write().

970  {
971  KALDI_ASSERT(v.Dim() == this->NumRows());
972  Real *data = this->data_;
973  const OtherReal *v_data = v.Data();
974  MatrixIndexT nr = this->num_rows_;
975  for (MatrixIndexT i = 0; i < nr; i++)
976  for (MatrixIndexT j = 0; j <= i; j++, data++)
977  *data += alpha * v_data[i] * v_data[j];
978 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void AddVec2 ( const double  alpha,
const VectorBase< double > &  v 
)
void AddVec2 ( const float  alpha,
const VectorBase< float > &  v 
)

Definition at line 939 of file sp-matrix.cc.

References kaldi::cblas_Xspr(), VectorBase< Real >::Data(), data_, VectorBase< Real >::Dim(), and KALDI_ASSERT.

939  {
940  KALDI_ASSERT(v.Dim() == this->NumRows());
941  cblas_Xspr(v.Dim(), alpha, v.Data(), 1,
942  this->data_);
943 }
MatrixIndexT NumRows() const
void cblas_Xspr(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, float *Adata)
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void AddVec2 ( const double  alpha,
const VectorBase< double > &  v 
)

Definition at line 962 of file sp-matrix.cc.

References kaldi::cblas_Xspr(), VectorBase< Real >::Data(), data_, VectorBase< Real >::Dim(), and KALDI_ASSERT.

962  {
963  KALDI_ASSERT(v.Dim() == num_rows_);
964  cblas_Xspr(v.Dim(), alpha, v.Data(), 1, data_);
965 }
void cblas_Xspr(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, float *Adata)
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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 at line 946 of file sp-matrix.cc.

References VectorBase< Real >::Data(), PackedMatrix< Real >::Data(), data_, VectorBase< Real >::Dim(), KALDI_ASSERT, and PackedMatrix< Real >::NumRows().

Referenced by kaldi::SolveQuadraticMatrixProblem(), kaldi::SolveQuadraticProblem(), and kaldi::UnitTestAddVec2Sp().

947  {
948  KALDI_ASSERT(v.Dim() == this->NumRows() && S.NumRows() == this->NumRows());
949  const Real *Sdata = S.Data();
950  const Real *vdata = v.Data();
951  Real *data = this->data_;
952  MatrixIndexT dim = this->num_rows_;
953  for (MatrixIndexT r = 0; r < dim; r++)
954  for (MatrixIndexT c = 0; c <= r; c++, Sdata++, data++)
955  *data = beta * *data + alpha * vdata[r] * vdata[c] * *Sdata;
956 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void AddVecVec ( const Real  alpha,
const VectorBase< Real > &  v,
const VectorBase< Real > &  w 
)

rank-two update, this <– this + alpha (v w' + w v').

Definition at line 1171 of file sp-matrix.cc.

References kaldi::cblas_Xspr2(), VectorBase< Real >::Data(), data_, VectorBase< Real >::Dim(), and KALDI_ASSERT.

Referenced by kaldi::UnitTestSpAddVecVec().

1172  {
1173  int32 dim = this->NumRows();
1174  KALDI_ASSERT(dim == v.Dim() && dim == w.Dim() && dim > 0);
1175  cblas_Xspr2(dim, alpha, v.Data(), 1, w.Data(), 1, this->data_);
1176 }
void cblas_Xspr2(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, const float *Ydata, MatrixIndexT incY, float *Adata)
MatrixIndexT NumRows() const
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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.

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.

Definition at line 560 of file sp-matrix.cc.

References SpMatrix< Real >::AddMat2(), SpMatrix< Real >::AddMat2Sp(), VectorBase< Real >::ApplyPow(), TpMatrix< Real >::Cholesky(), VectorBase< Real >::Dim(), SpMatrix< Real >::Eig(), rnnlm::i, TpMatrix< Real >::Invert(), KALDI_ASSERT, KALDI_LOG, kaldi::kNoTrans, MatrixBase< Real >::MulColsVec(), PackedMatrix< Real >::NumRows(), and PackedMatrix< Real >::Scale().

Referenced by main(), kaldi::UnitTestFloorChol(), kaldi::UnitTestFloorUnit(), IvectorExtractorStats::UpdateVariances(), and EbwAmSgmm2Updater::UpdateVars().

561  {
562  MatrixIndexT dim = this->NumRows();
563  int nfloored = 0;
564  KALDI_ASSERT(C.NumRows() == dim);
565  KALDI_ASSERT(alpha > 0);
566  TpMatrix<Real> L(dim);
567  L.Cholesky(C);
568  L.Scale(std::sqrt(alpha)); // equivalent to scaling C by alpha.
569  TpMatrix<Real> LInv(L);
570  LInv.Invert();
571 
572  SpMatrix<Real> D(dim);
573  { // D = L^{-1} * (*this) * L^{-T}
574  Matrix<Real> LInvFull(LInv);
575  D.AddMat2Sp(1.0, LInvFull, kNoTrans, (*this), 0.0);
576  }
577 
578  Vector<Real> l(dim);
579  Matrix<Real> U(dim, dim);
580 
581  D.Eig(&l, &U);
582 
583  if (verbose) {
584  KALDI_LOG << "ApplyFloor: flooring following diagonal to 1: " << l;
585  }
586  for (MatrixIndexT i = 0; i < l.Dim(); i++) {
587  if (l(i) < 1.0) {
588  nfloored++;
589  l(i) = 1.0;
590  }
591  }
592  l.ApplyPow(0.5);
593  U.MulColsVec(l);
594  D.AddMat2(1.0, U, kNoTrans, 0.0);
595  { // D' := U * diag(l') * U^T ... l'=floor(l, 1)
596  Matrix<Real> LFull(L);
597  (*this).AddMat2Sp(1.0, LFull, kNoTrans, D, 0.0); // A := L * D' * L^T
598  }
599  return nfloored;
600 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
#define KALDI_LOG
Definition: kaldi-error.h:133
int ApplyFloor ( Real  floor)

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.

Definition at line 613 of file sp-matrix.cc.

References rnnlm::i, and kaldi::kNoTrans.

613  {
614  MatrixIndexT Dim = this->NumRows();
615  int nfloored = 0;
616  Vector<Real> s(Dim);
617  Matrix<Real> P(Dim, Dim);
618  (*this).Eig(&s, &P);
619  for (MatrixIndexT i = 0; i < Dim; i++) {
620  if (s(i) < floor) {
621  nfloored++;
622  s(i) = floor;
623  }
624  }
625  (*this).AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
626  return nfloored;
627 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
void ApplyPow ( Real  exponent)

Takes matrix to a fraction power via Svd.

Will throw exception if matrix is not positive semidefinite (to within a tolerance)

Definition at line 117 of file sp-matrix.cc.

References VectorBase< Real >::ApplyPow(), KALDI_ASSERT, KALDI_ERR, kaldi::kNoTrans, and MatrixBase< Real >::MulColsVec().

Referenced by AmSgmm2::SplitSubstates(), and kaldi::UnitTestPower().

117  {
118  if (power == 1) return; // can do nothing.
119  MatrixIndexT D = this->NumRows();
120  KALDI_ASSERT(D > 0);
121  Matrix<Real> U(D, D);
122  Vector<Real> l(D);
123  (*this).SymPosSemiDefEig(&l, &U);
124 
125  Vector<Real> l_copy(l);
126  try {
127  l.ApplyPow(power * 0.5);
128  }
129  catch(...) {
130  KALDI_ERR << "Error taking power " << (power * 0.5) << " of vector "
131  << l_copy;
132  }
133  U.MulColsVec(l);
134  (*this).AddMat2(1.0, U, kNoTrans, 0.0);
135 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ERR
Definition: kaldi-error.h:127
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
bool ApproxEqual ( const SpMatrix< Real > &  other,
float  tol = 0.01 
) const

Returns true if ((*this)-other).FrobeniusNorm() <= tol*(*this).FrobeniusNorma()

Definition at line 549 of file sp-matrix.cc.

References SpMatrix< Real >::AddSp(), SpMatrix< Real >::FrobeniusNorm(), KALDI_ERR, and PackedMatrix< Real >::NumRows().

Referenced by kaldi::UnitTestEigSp(), and kaldi::UnitTestNorm().

549  {
550  if (this->NumRows() != other.NumRows())
551  KALDI_ERR << "SpMatrix::AproxEqual, size mismatch, "
552  << this->NumRows() << " vs. " << other.NumRows();
553  SpMatrix<Real> tmp(*this);
554  tmp.AddSp(-1.0, other);
555  return (tmp.FrobeniusNorm() <= tol * std::max(this->FrobeniusNorm(), other.FrobeniusNorm()));
556 }
Real FrobeniusNorm() const
sqrt of sum of square elements.
Definition: sp-matrix.cc:537
MatrixIndexT NumRows() const
#define KALDI_ERR
Definition: kaldi-error.h:127
Real Cond ( ) const
inline

Returns maximum ratio of singular values.

Definition at line 145 of file sp-matrix.h.

Referenced by IvectorExtractor::GetIvectorDistribution(), kaldi::InitRand(), and kaldi::InitRandNonsingular().

145  {
146  Matrix<Real> tmp(*this);
147  return tmp.Cond();
148  }
void CopyFromMat ( const MatrixBase< Real > &  orig,
SpCopyType  copy_type = kTakeMean 
)

Definition at line 138 of file sp-matrix.cc.

References MatrixBase< Real >::Data(), data_, rnnlm::i, rnnlm::j, KALDI_ASSERT, KALDI_ERR, kaldi::kTakeLower, kaldi::kTakeMean, kaldi::kTakeMeanAndCheck, kaldi::kTakeUpper, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), and MatrixBase< Real >::Stride().

Referenced by MatrixExponential< Real >::Backprop(), MatrixExponential< Real >::Compute(), BasisFmllrEstimate::ComputeAmDiagPrecond(), SpMatrix< BaseFloat >::SpMatrix(), kaldi::UnitTestCopySp(), kaldi::UnitTestDeterminant(), kaldi::UnitTestDeterminantSign(), and kaldi::UnitTestPower().

139  {
140  KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
141  MatrixIndexT D = this->NumRows();
142 
143  switch (copy_type) {
144  case kTakeMeanAndCheck:
145  {
146  Real good_sum = 0.0, bad_sum = 0.0;
147  for (MatrixIndexT i = 0; i < D; i++) {
148  for (MatrixIndexT j = 0; j < i; j++) {
149  Real a = M(i, j), b = M(j, i), avg = 0.5*(a+b), diff = 0.5*(a-b);
150  (*this)(i, j) = avg;
151  good_sum += std::abs(avg);
152  bad_sum += std::abs(diff);
153  }
154  good_sum += std::abs(M(i, i));
155  (*this)(i, i) = M(i, i);
156  }
157  if (bad_sum > 0.01 * good_sum) {
158  KALDI_ERR << "SpMatrix::Copy(), source matrix is not symmetric: "
159  << bad_sum << ">" << good_sum;
160  }
161  break;
162  }
163  case kTakeMean:
164  {
165  for (MatrixIndexT i = 0; i < D; i++) {
166  for (MatrixIndexT j = 0; j < i; j++) {
167  (*this)(i, j) = 0.5*(M(i, j) + M(j, i));
168  }
169  (*this)(i, i) = M(i, i);
170  }
171  break;
172  }
173  case kTakeLower:
174  { // making this one a bit more efficient.
175  const Real *src = M.Data();
176  Real *dest = this->data_;
177  MatrixIndexT stride = M.Stride();
178  for (MatrixIndexT i = 0; i < D; i++) {
179  for (MatrixIndexT j = 0; j <= i; j++)
180  dest[j] = src[j];
181  dest += i + 1;
182  src += stride;
183  }
184  }
185  break;
186  case kTakeUpper:
187  for (MatrixIndexT i = 0; i < D; i++)
188  for (MatrixIndexT j = 0; j <= i; j++)
189  (*this)(i, j) = M(j, i);
190  break;
191  default:
192  KALDI_ASSERT("Invalid argument to SpMatrix::CopyFromMat");
193  }
194 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ERR
Definition: kaldi-error.h:127
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void CopyFromSp ( const SpMatrix< OtherReal > &  other)
inline

Definition at line 90 of file sp-matrix.h.

90  {
92  }
void CopyFromPacked(const PackedMatrix< OtherReal > &orig)
template 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.

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.

Definition at line 433 of file qr.cc.

References VectorBase< Real >::CopyDiagFromPacked(), VectorBase< Real >::Dim(), KALDI_ASSERT, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), SpMatrix< Real >::Qr(), MatrixBase< Real >::Transpose(), and SpMatrix< Real >::Tridiagonalize().

Referenced by SpMatrix< Real >::ApplyFloor(), kaldi::ComputeFeatureNormalizingTransform(), AmSgmm2::ComputeFmllrPreXform(), kaldi::ComputeLdaTransform(), kaldi::ComputePca(), kaldi::EstimateSgmm2FmllrSubspace(), kaldi::GetLogDetNoFailure(), PldaEstimator::GetOutput(), IvectorExtractor::InvertWithFlooring(), main(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), OnlineNaturalGradient::PreconditionDirectionsInternal(), OnlinePreconditioner::PreconditionDirectionsInternal(), kaldi::UnitTestEigSp(), kaldi::UnitTestMaxAbsEig(), kaldi::UnitTestPca2(), kaldi::UnitTestPldaEstimation(), kaldi::UnitTestSvdSpeed(), kaldi::UnitTestTopEigs(), PldaUnsupervisedAdaptor::UpdatePlda(), and IvectorExtractorStats::UpdatePrior().

433  {
434  MatrixIndexT dim = this->NumRows();
435  KALDI_ASSERT(s->Dim() == dim);
436  KALDI_ASSERT(P == NULL || (P->NumRows() == dim && P->NumCols() == dim));
437 
438  SpMatrix<Real> A(*this); // Copy *this, since the tridiagonalization
439  // and QR decomposition are destructive.
440  // Note: for efficiency of memory access, the tridiagonalization
441  // algorithm makes the *rows* of P the eigenvectors, not the columns.
442  // We'll transpose P before we exit.
443  // Also note: P may be null if you don't want the eigenvectors. This
444  // will make this function more efficient.
445 
446  A.Tridiagonalize(P); // Tridiagonalizes.
447  A.Qr(P); // Diagonalizes.
448  if(P) P->Transpose();
449  s->CopyDiagFromPacked(A);
450 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void EigInternal ( VectorBase< Real > *  s,
MatrixBase< Real > *  P,
Real  tolerance,
int  recurse 
) const
private
void Exp ( )

Definition at line 85 of file sp-matrix.cc.

References MatrixExponential< Real >::Compute(), and KALDI_ASSERT.

Referenced by kaldi::UnitTestSpLogExp().

85  {
86  // The most natural way to do this would be to do a symmetric eigenvalue
87  // decomposition, but in order to work with basic ATLAS without CLAPACK, we
88  // don't have symmetric eigenvalue decomposition code (correction: didn't have
89  // it). Instead we use the MatrixExponential class which expands it out as a
90  // Taylor series (with a pre-scaling trick to make it reasonably fast).
91 
92  KALDI_ASSERT(this->NumRows() != 0);
93  Matrix<Real> M(*this), expM(this->NumRows(), this->NumRows());
94  MatrixExponential<Real> me;
95  me.Compute(M, &expM); // compute exp(M)
96  this->CopyFromMat(expM); // by default, checks it's symmetric.
97 }
MatrixIndexT NumRows() const
void CopyFromMat(const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
Definition: sp-matrix.cc:138
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
Real FrobeniusNorm ( ) const

sqrt of sum of square elements.

Definition at line 537 of file sp-matrix.cc.

References rnnlm::i, and rnnlm::j.

Referenced by SpMatrix< Real >::ApproxEqual(), kaldi::UnitTestCuSpMatrixApproxEqual(), kaldi::UnitTestCuSpMatrixSetUnit(), and kaldi::UnitTestNorm().

537  {
538  Real sum = 0.0;
539  MatrixIndexT R = this->NumRows();
540  for (MatrixIndexT i = 0; i < R; i++) {
541  for (MatrixIndexT j = 0; j < i; j++)
542  sum += (*this)(i, j) * (*this)(i, j) * 2;
543  sum += (*this)(i, i) * (*this)(i, i);
544  }
545  return std::sqrt(sum);
546 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
void Invert ( Real *  logdet = NULL,
Real *  det_sign = NULL,
bool  inverse_needed = true 
)

matrix inverse.

if inverse_needed = false, will fill matrix with garbage. (only useful if logdet wanted).

Definition at line 245 of file sp-matrix.cc.

References kaldi::clapack_Xsptrf(), kaldi::clapack_Xsptri(), data_, kaldi::diag, rnnlm::i, KALDI_ASSERT, KALDI_ERR, KALDI_MEMALIGN, KALDI_MEMALIGN_FREE, and kaldi::Log().

Referenced by kaldi::ComputeMllrMatrix(), PldaEstimator::ComputeObjfPart1(), PldaEstimator::ComputeObjfPart2(), DiagGmm::CopyFromFullGmm(), IvectorExtractor::GetIvectorDistribution(), GetLogLikeTest(), PldaEstimator::GetStatsFromClassMeans(), SpMatrix< Real >::InvertDouble(), SpMatrix< Real >::LogDet(), main(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), Sgmm2Project::ProjectVariance(), kaldi::UnitTestCuMatrixSymInvertPosDef(), kaldi::UnitTestCuSpMatrixInvert(), kaldi::UnitTestLbfgs(), kaldi::UnitTestSolve(), kaldi::UnitTestSpInvert(), IvectorExtractorStats::UpdateVariances(), and EbwAmSgmm2Updater::UpdateVars().

245  {
246  // these are CLAPACK types
247  KaldiBlasInt result;
248  KaldiBlasInt rows = static_cast<int>(this->num_rows_);
249  KaldiBlasInt* p_ipiv = new KaldiBlasInt[rows];
250  Real *p_work; // workspace for the lapack function
251  void *temp;
252  if ((p_work = static_cast<Real*>(
253  KALDI_MEMALIGN(16, sizeof(Real) * rows, &temp))) == NULL)
254  throw std::bad_alloc();
255 #ifdef HAVE_OPENBLAS
256  memset(p_work, 0, sizeof(Real) * rows); // gets rid of a probably
257  // spurious Valgrind warning about jumps depending upon uninitialized values.
258 #endif
259 
260 
261  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
262  // of the data. We have a row-wise storage, therefore, we need to "invert"
263  clapack_Xsptrf(&rows, this->data_, p_ipiv, &result);
264 
265 
266  KALDI_ASSERT(result >= 0 && "Call to CLAPACK ssptrf_ called with wrong arguments");
267 
268  if (result > 0) { // Singular...
269  if (det_sign) *det_sign = 0;
270  if (logdet) *logdet = -std::numeric_limits<Real>::infinity();
271  if (need_inverse) KALDI_ERR << "CLAPACK stptrf_ : factorization failed";
272  } else { // Not singular.. compute log-determinant if needed.
273  if (logdet != NULL || det_sign != NULL) {
274  Real prod = 1.0, log_prod = 0.0;
275  int sign = 1;
276  for (int i = 0; i < (int)this->num_rows_; i++) {
277  if (p_ipiv[i] > 0) { // not a 2x2 block...
278  // if (p_ipiv[i] != i+1) sign *= -1; // row swap.
279  Real diag = (*this)(i, i);
280  prod *= diag;
281  } else { // negative: 2x2 block. [we are in first of the two].
282  i++; // skip over the first of the pair.
283  // each 2x2 block...
284  Real diag1 = (*this)(i, i), diag2 = (*this)(i-1, i-1),
285  offdiag = (*this)(i, i-1);
286  Real thisdet = diag1*diag2 - offdiag*offdiag;
287  // thisdet == determinant of 2x2 block.
288  // The following line is more complex than it looks: there are 2 offsets of
289  // 1 that cancel.
290  prod *= thisdet;
291  }
292  if (i == (int)(this->num_rows_-1) || fabs(prod) < 1.0e-10 || fabs(prod) > 1.0e+10) {
293  if (prod < 0) { prod = -prod; sign *= -1; }
294  log_prod += kaldi::Log(std::abs(prod));
295  prod = 1.0;
296  }
297  }
298  if (logdet != NULL) *logdet = log_prod;
299  if (det_sign != NULL) *det_sign = sign;
300  }
301  }
302  if (!need_inverse) {
303  delete [] p_ipiv;
304  KALDI_MEMALIGN_FREE(p_work);
305  return; // Don't need what is computed next.
306  }
307  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
308  // of the data. We have a row-wise storage, therefore, we need to "invert"
309  clapack_Xsptri(&rows, this->data_, p_ipiv, p_work, &result);
310 
311  KALDI_ASSERT(result >=0 &&
312  "Call to CLAPACK ssptri_ called with wrong arguments");
313 
314  if (result != 0) {
315  KALDI_ERR << "CLAPACK ssptrf_ : Matrix is singular";
316  }
317 
318  delete [] p_ipiv;
319  KALDI_MEMALIGN_FREE(p_work);
320 }
void clapack_Xsptri(KaldiBlasInt *num_rows, float *Mdata, KaldiBlasInt *ipiv, float *work, KaldiBlasInt *result)
double Log(double x)
Definition: kaldi-math.h:100
#define KALDI_MEMALIGN(align, size, pp_orig)
Definition: kaldi-utils.h:54
#define KALDI_ERR
Definition: kaldi-error.h:127
#define KALDI_MEMALIGN_FREE(x)
Definition: kaldi-utils.h:56
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void clapack_Xsptrf(KaldiBlasInt *num_rows, float *Mdata, KaldiBlasInt *ipiv, KaldiBlasInt *result)
void InvertDouble ( Real *  logdet = NULL,
Real *  det_sign = NULL,
bool  inverse_needed = true 
)

Definition at line 336 of file sp-matrix.cc.

References SpMatrix< Real >::Invert().

Referenced by AmSgmm2::ComputeFmllrPreXform(), FullGmm::ComputeGconsts(), FullGmm::GetComponentMean(), FullGmm::GetMeans(), main(), FullGmm::SetInvCovars(), and kaldi::UnitTestSpInvert().

337  {
338  SpMatrix<double> dmat(*this);
339  double logdet_tmp, det_sign_tmp;
340  dmat.Invert(logdet ? &logdet_tmp : NULL,
341  det_sign ? &det_sign_tmp : NULL,
342  inverse_needed);
343  if (logdet) *logdet = logdet_tmp;
344  if (det_sign) *det_sign = det_sign_tmp;
345  (*this).CopyFromSp(dmat);
346 }
bool IsDiagonal ( Real  cutoff = 1.0e-05) const

Definition at line 489 of file sp-matrix.cc.

References rnnlm::i, and rnnlm::j.

Referenced by AmSgmm2::ComputeFmllrPreXform(), MleAmSgmm2Updater::RenormalizeV(), kaldi::SolveDoubleQuadraticMatrixProblem(), kaldi::UnitTestPca2(), and kaldi::UnitTestTridiagonalizeAndQr().

489  {
490  MatrixIndexT R = this->NumRows();
491  Real bad_sum = 0.0, good_sum = 0.0;
492  for (MatrixIndexT i = 0; i < R; i++) {
493  for (MatrixIndexT j = 0; j <= i; j++) {
494  if (i == j)
495  good_sum += std::abs((*this)(i, j));
496  else
497  bad_sum += std::abs((*this)(i, j));
498  }
499  }
500  return (!(bad_sum > good_sum * cutoff));
501 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
bool IsPosDef ( ) const

Definition at line 101 of file sp-matrix.cc.

References TpMatrix< Real >::Cholesky(), and KALDI_ASSERT.

Referenced by MleAmSgmm2Updater::RenormalizeV().

101  {
102  MatrixIndexT D = (*this).NumRows();
103  KALDI_ASSERT(D > 0);
104  try {
105  TpMatrix<Real> C(D);
106  C.Cholesky(*this);
107  for (MatrixIndexT r = 0; r < D; r++)
108  if (C(r, r) == 0.0) return false;
109  return true;
110  }
111  catch(...) { // not positive semidefinite.
112  return false;
113  }
114 }
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
bool IsTridiagonal ( Real  cutoff = 1.0e-05) const

Definition at line 515 of file sp-matrix.cc.

References rnnlm::i, and rnnlm::j.

Referenced by kaldi::UnitTestTridiag().

515  {
516  MatrixIndexT R = this->NumRows();
517  Real max_abs_2diag = 0.0, max_abs_offdiag = 0.0;
518  for (MatrixIndexT i = 0; i < R; i++)
519  for (MatrixIndexT j = 0; j <= i; j++) {
520  if (j+1 < i)
521  max_abs_offdiag = std::max(max_abs_offdiag,
522  std::abs((*this)(i, j)));
523  else
524  max_abs_2diag = std::max(max_abs_2diag,
525  std::abs((*this)(i, j)));
526  }
527  return (max_abs_offdiag <= cutoff * max_abs_2diag);
528 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
bool IsZero ( Real  cutoff = 1.0e-05) const

Definition at line 531 of file sp-matrix.cc.

Referenced by kaldi::GpsrBasic(), kaldi::GpsrBB(), kaldi::SolveQuadraticMatrixProblem(), kaldi::SolveQuadraticProblem(), and kaldi::UnitTestResize().

531  {
532  if (this->num_rows_ == 0) return true;
533  return (this->Max() <= cutoff && this->Min() >= -cutoff);
534 }
Real Min() const
MatrixIndexT num_rows_
Real Max() const
MatrixIndexT LimitCond ( Real  maxCond = 1.0e+5,
bool  invert = false 
)

Definition at line 630 of file sp-matrix.cc.

References rnnlm::i, KALDI_ASSERT, KALDI_WARN, kaldi::kNoTrans, VectorBase< Real >::Max(), and MatrixBase< Real >::MulColsVec().

Referenced by MleAmSgmm2Updater::ComputeMPrior(), SpMatrix< BaseFloat >::LimitCondDouble(), kaldi::UnitTestLimitCond(), and kaldi::UnitTestLimitCondInvert().

630  { // e.g. maxCond = 1.0e+05.
631  MatrixIndexT Dim = this->NumRows();
632  Vector<Real> s(Dim);
633  Matrix<Real> P(Dim, Dim);
634  (*this).SymPosSemiDefEig(&s, &P);
635  KALDI_ASSERT(maxCond > 1);
636  Real floor = s.Max() / maxCond;
637  if (floor < 0) floor = 0;
638  if (floor < 1.0e-40) {
639  KALDI_WARN << "LimitCond: limiting " << floor << " to 1.0e-40";
640  floor = 1.0e-40;
641  }
642  MatrixIndexT nfloored = 0;
643  for (MatrixIndexT i = 0; i < Dim; i++) {
644  if (s(i) <= floor) nfloored++;
645  if (invert)
646  s(i) = 1.0 / std::sqrt(std::max(s(i), floor));
647  else
648  s(i) = std::sqrt(std::max(s(i), floor));
649  }
650  P.MulColsVec(s);
651  (*this).AddMat2(1.0, P, kNoTrans, 0.0); // (*this) = P*P^T. ... (*this) = P * floor(s) * P^T ... if P was original P.
652  return nfloored;
653 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_WARN
Definition: kaldi-error.h:130
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
MatrixIndexT LimitCondDouble ( Real  maxCond = 1.0e+5,
bool  invert = false 
)
inline

Definition at line 343 of file sp-matrix.h.

Referenced by AmSgmm2::ComputeHsmFromModel(), and MleAmSgmm2Updater::UpdateVars().

343  {
344  SpMatrix<double> dmat(*this);
345  MatrixIndexT ans = dmat.LimitCond(maxCond, invert);
346  (*this).CopyFromSp(dmat);
347  return ans;
348  }
int32 MatrixIndexT
Definition: matrix-common.h:96
void Log ( )

Takes log of the matrix (does eigenvalue decomposition then takes log of eigenvalues and reconstructs).

Will throw of not +ve definite.

Definition at line 74 of file sp-matrix.cc.

References VectorBase< Real >::ApplyLog(), KALDI_ASSERT, and kaldi::kNoTrans.

Referenced by kaldi::UnitTestSpLogExp().

74  {
75  KALDI_ASSERT(this->NumRows() != 0);
76  Vector<Real> s(this->NumRows());
77  Matrix<Real> P(this->NumRows(), this->NumRows());
78  SymPosSemiDefEig(&s, &P);
79  s.ApplyLog(); // Per-element log.
80  // this <-- P * diag(s) * P^T
81  this->AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
82 }
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:1105
MatrixIndexT NumRows() const
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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
Real LogDet ( Real *  det_sign = NULL) const

Definition at line 603 of file sp-matrix.cc.

References SpMatrix< Real >::Invert().

Referenced by kaldi::UnitTestDeterminant(), kaldi::UnitTestDeterminantSign(), kaldi::UnitTestTridiagonalize(), kaldi::UnitTestTridiagonalizeAndQr(), and EbwAmSgmm2Updater::UpdateVars().

603  {
604  Real log_det;
605  SpMatrix<Real> tmp(*this);
606  // false== output not needed (saves some computation).
607  tmp.Invert(&log_det, det_sign, false);
608  return log_det;
609 }
Real LogPosDefDet ( ) const

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()

Definition at line 36 of file sp-matrix.cc.

References TpMatrix< Real >::Cholesky(), kaldi::diag, rnnlm::i, and kaldi::Log().

Referenced by FullGmm::ComputeGconsts(), MleAmSgmm2Updater::ComputeMPrior(), kaldi::GetLogDetNoFailure(), FullGmm::MergedComponentsLogdet(), IvectorExtractorStats::PriorDiagnostics(), rand_posdef_spmatrix(), kaldi::unittest::RandPosdefSpMatrix(), RandPosdefSpMatrix(), kaldi::UnitTestDeterminant(), kaldi::UnitTestDeterminantSign(), and IvectorExtractorStats::UpdateVariances().

36  {
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 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
double Log(double x)
Definition: kaldi-math.h:100
Real MaxAbsEig ( ) const

Returns the maximum of the absolute values of any of the eigenvalues.

Definition at line 67 of file sp-matrix.cc.

References VectorBase< Real >::Max(), and VectorBase< Real >::Min().

Referenced by main(), kaldi::UnitTestLinearCgd(), kaldi::UnitTestMaxAbsEig(), and IvectorExtractorStats::UpdateVariances().

67  {
68  Vector<Real> s(this->NumRows());
69  this->Eig(&s, static_cast<MatrixBase<Real>*>(NULL));
70  return std::max(s.Max(), -s.Min());
71 }
MatrixIndexT NumRows() const
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
Real operator() ( MatrixIndexT  r,
MatrixIndexT  c 
) const
inline

Definition at line 102 of file sp-matrix.h.

102  {
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  }
void swap(basic_filebuf< CharT, Traits > &x, basic_filebuf< CharT, Traits > &y)
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
Real& operator() ( MatrixIndexT  r,
MatrixIndexT  c 
)
inline

Definition at line 115 of file sp-matrix.h.

115  {
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  }
void swap(basic_filebuf< CharT, Traits > &x, basic_filebuf< CharT, Traits > &y)
MatrixIndexT num_rows_
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
SpMatrix<Real>& operator= ( const SpMatrix< Real > &  other)
inline

Definition at line 126 of file sp-matrix.h.

126  {
128  return *this;
129  }
PackedMatrix< Real > & operator=(const PackedMatrix< Real > &other)
Definition: packed-matrix.h:68
void PrintEigs ( const char *  name)
inline

Definition at line 213 of file sp-matrix.h.

Referenced by MleAmSgmm2Updater::ComputeMPrior(), kaldi::nnet2::PreconditionDirections(), and MleAmSgmm2Updater::RenormalizeV().

213  {
214  Vector<Real> s((*this).NumRows());
215  Matrix<Real> P((*this).NumRows(), (*this).NumCols());
216  SymPosSemiDefEig(&s, &P);
217  KALDI_LOG << "PrintEigs: " << name << ": " << s;
218  }
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
#define KALDI_LOG
Definition: kaldi-error.h:133
template void Qr ( MatrixBase< Real > *  Q)

The symmetric QR algorithm.

This is the symmetric QR algorithm, from Golub and Van Loan 3rd ed., Algorithm 8.3.3.

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.

Q is transposed w.r.t. there, though.

Definition at line 409 of file qr.cc.

References VectorBase< Real >::Data(), kaldi::diag, rnnlm::i, KALDI_ASSERT, rnnlm::n, MatrixBase< Real >::NumRows(), and kaldi::QrInternal().

Referenced by SpMatrix< Real >::Eig(), SpMatrix< Real >::TopEigs(), and kaldi::UnitTestTridiagonalizeAndQr().

409  {
410  KALDI_ASSERT(this->IsTridiagonal());
411  // We envisage that Q would be square but we don't check for this,
412  // as there are situations where you might not want this.
413  KALDI_ASSERT(Q == NULL || Q->NumRows() == this->NumRows());
414  // Note: the first couple of lines of the algorithm they give would be done
415  // outside of this function, by calling Tridiagonalize().
416 
417  MatrixIndexT n = this->NumRows();
418  Vector<Real> diag(n), off_diag(n-1);
419  for (MatrixIndexT i = 0; i < n; i++) {
420  diag(i) = (*this)(i, i);
421  if (i > 0) off_diag(i-1) = (*this)(i, i-1);
422  }
423  QrInternal(n, diag.Data(), off_diag.Data(), Q);
424  // Now set *this to the value represented by diag and off_diag.
425  this->SetZero();
426  for (MatrixIndexT i = 0; i < n; i++) {
427  (*this)(i, i) = diag(i);
428  if (i > 0) (*this)(i, i-1) = off_diag(i-1);
429  }
430 }
void QrInternal(MatrixIndexT n, Real *diag, Real *off_diag, MatrixBase< Real > *Q)
Definition: qr.cc:331
bool IsTridiagonal(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:515
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
struct rnnlm::@11::@12 n
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void Swap ( SpMatrix< Real > *  other)

Shallow swap.

Definition at line 51 of file sp-matrix.cc.

References data_, PackedMatrix< Real >::data_, PackedMatrix< Real >::num_rows_, and kaldi::swap().

51  {
52  std::swap(this->data_, other->data_);
53  std::swap(this->num_rows_, other->num_rows_);
54 }
void swap(basic_filebuf< CharT, Traits > &x, basic_filebuf< CharT, Traits > &y)
MatrixIndexT num_rows_
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.

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.

Definition at line 57 of file sp-matrix.cc.

References VectorBase< Real >::ApplyFloor(), KALDI_ASSERT, VectorBase< Real >::Max(), and VectorBase< Real >::Min().

Referenced by Sgmm2Project::ComputeLdaTransform(), BasisFmllrEstimate::EstimateFmllrBasis(), GetLogLikeTest(), SpMatrix< BaseFloat >::PrintEigs(), MleAmSgmm2Updater::RenormalizeN(), MleAmSgmm2Updater::RenormalizeV(), kaldi::SolveDoubleQuadraticMatrixProblem(), kaldi::SolveQuadraticMatrixProblem(), kaldi::SolveQuadraticProblem(), kaldi::UnitTestFloorUnit(), kaldi::UnitTestMat2Vec(), and kaldi::UnitTestSvdBad().

59  {
60  Eig(s, P);
61  Real max = s->Max(), min = s->Min();
62  KALDI_ASSERT(-min <= tolerance * max);
63  s->ApplyFloor(0.0);
64 }
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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
template 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 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.

Definition at line 454 of file qr.cc.

References MatrixBase< Real >::AddMatMat(), VectorBase< Real >::AddSpVec(), VectorBase< Real >::AddVec(), VectorBase< Real >::CopyDiagFromSp(), MatrixBase< Real >::CopyFromMat(), VectorBase< Real >::CopyFromVec(), rnnlm::d, VectorBase< Real >::Dim(), KALDI_ASSERT, KALDI_ERR, kaldi::kTrans, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), SpMatrix< Real >::Qr(), VectorBase< Real >::Range(), MatrixBase< Real >::Range(), VectorBase< Real >::Scale(), VectorBase< Real >::SetRandn(), MatrixBase< Real >::SetUnit(), kaldi::SortSvd(), and kaldi::VecVec().

Referenced by kaldi::ComputePca(), and kaldi::UnitTestTopEigs().

455  {
456  const SpMatrix<Real> &S(*this); // call this "S" for easy notation.
457  MatrixIndexT eig_dim = s->Dim(); // Space of dim we want to retain.
458  if (lanczos_dim <= 0)
459  lanczos_dim = std::max(eig_dim + 50, eig_dim + eig_dim/2);
460  MatrixIndexT dim = this->NumRows();
461  if (lanczos_dim >= dim) {
462  // There would be no speed advantage in using this method, so just
463  // use the regular approach.
464  Vector<Real> s_tmp(dim);
465  Matrix<Real> P_tmp(dim, dim);
466  this->Eig(&s_tmp, &P_tmp);
467  SortSvd(&s_tmp, &P_tmp);
468  s->CopyFromVec(s_tmp.Range(0, eig_dim));
469  P->CopyFromMat(P_tmp.Range(0, dim, 0, eig_dim));
470  return;
471  }
472  KALDI_ASSERT(eig_dim <= dim && eig_dim > 0);
473  KALDI_ASSERT(P->NumRows() == dim && P->NumCols() == eig_dim); // each column
474  // is one eigenvector.
475 
476  Matrix<Real> Q(lanczos_dim, dim); // The rows of Q will be the
477  // orthogonal vectors of the Krylov subspace.
478 
479  SpMatrix<Real> T(lanczos_dim); // This will be equal to Q S Q^T,
480  // i.e. *this projected into the Krylov subspace. Note: only the
481  // diagonal and off-diagonal fo T are nonzero, i.e. it's tridiagonal,
482  // but we don't have access to the low-level algorithms that work
483  // on that type of matrix (since we want to use ATLAS). So we just
484  // do normal SVD, on a full matrix; it won't typically dominate.
485 
486  Q.Row(0).SetRandn();
487  Q.Row(0).Scale(1.0 / Q.Row(0).Norm(2));
488  for (MatrixIndexT d = 0; d < lanczos_dim; d++) {
489  Vector<Real> r(dim);
490  r.AddSpVec(1.0, S, Q.Row(d), 0.0);
491  // r = S * q_d
492  MatrixIndexT counter = 0;
493  Real end_prod;
494  while (1) { // Normally we'll do this loop only once:
495  // we repeat to handle cases where r gets very much smaller
496  // and we want to orthogonalize again.
497  // We do "full orthogonalization" to preserve stability,
498  // even though this is usually a waste of time.
499  Real start_prod = VecVec(r, r);
500  for (SignedMatrixIndexT e = d; e >= 0; e--) { // e must be signed!
501  SubVector<Real> q_e(Q, e);
502  Real prod = VecVec(r, q_e);
503  if (counter == 0 && static_cast<MatrixIndexT>(e) + 1 >= d) // Keep T tridiagonal, which
504  T(d, e) = prod; // mathematically speaking, it is.
505  r.AddVec(-prod, q_e); // Subtract component in q_e.
506  }
507  if (d+1 == lanczos_dim) break;
508  end_prod = VecVec(r, r);
509  if (end_prod <= 0.1 * start_prod) {
510  // also handles case where both are 0.
511  // We're not confident any more that it's completely
512  // orthogonal to the rest so we want to re-do.
513  if (end_prod == 0.0)
514  r.SetRandn(); // "Restarting".
515  counter++;
516  if (counter > 100)
517  KALDI_ERR << "Loop detected in Lanczos iteration.";
518  } else {
519  break;
520  }
521  }
522  if (d+1 != lanczos_dim) {
523  // OK, at this point we're satisfied that r is orthogonal
524  // to all previous rows.
525  KALDI_ASSERT(end_prod != 0.0); // should have looped.
526  r.Scale(1.0 / std::sqrt(end_prod)); // make it unit.
527  Q.Row(d+1).CopyFromVec(r);
528  }
529  }
530 
531  Matrix<Real> R(lanczos_dim, lanczos_dim);
532  R.SetUnit();
533  T.Qr(&R); // Diagonalizes T.
534  Vector<Real> s_tmp(lanczos_dim);
535  s_tmp.CopyDiagFromSp(T);
536 
537  // Now T = R * diag(s_tmp) * R^T.
538  // The next call sorts the elements of s from greatest to least absolute value,
539  // and moves around the rows of R in the corresponding way. This picks out
540  // the largest (absolute) eigenvalues.
541  SortSvd(&s_tmp, static_cast<Matrix<Real>*>(NULL), &R);
542  // Keep only the initial rows of R, those corresponding to greatest (absolute)
543  // eigenvalues.
544  SubMatrix<Real> Rsub(R, 0, eig_dim, 0, lanczos_dim);
545  SubVector<Real> s_sub(s_tmp, 0, eig_dim);
546  s->CopyFromVec(s_sub);
547 
548  // For working out what to do now, just assume the other eigenvalues were
549  // zero. This is just for purposes of knowing how to get the result, and
550  // not getting things wrongly transposed.
551  // We have T = Rsub^T * diag(s_sub) * Rsub.
552  // Now, T = Q S Q^T, with Q orthogonal, so S = Q^T T Q = Q^T Rsub^T * diag(s) * Rsub * Q.
553  // The output is P and we want S = P * diag(s) * P^T, so we need P = Q^T Rsub^T.
554  P->AddMatMat(1.0, Q, kTrans, Rsub, kTrans, 0.0);
555 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ERR
Definition: kaldi-error.h:127
int32 SignedMatrixIndexT
Definition: matrix-common.h:97
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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
Real VecVec(const VectorBase< Real > &a, const VectorBase< Real > &b)
Returns dot product between v1 and v2.
Definition: kaldi-vector.cc:36
void SortSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt, bool sort_on_absolute_value)
Function to ensure that SVD is sorted.
template void Tridiagonalize ( MatrixBase< Real > *  Q)

Tridiagonalize the matrix with an orthogonal transformation.

This routine tridiagonalizes *this.

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).

C.f. Golub and Van Loan 3rd ed., sec. 8.3.1 (p415). We reverse the order of the indices as it's more natural with packed lower-triangular matrices to do it this way. There's also a shift from one-based to zero-based indexing, so the index k is transformed k -> n - k, and a corresponding transpose...

Let the original *this be A. This algorithms replaces *this with a tridiagonal matrix T such that T = Q A Q^T for an orthogonal Q. Caution: Q is transposed vs. Golub and Van Loan. If Q != NULL it outputs Q.

Definition at line 165 of file qr.cc.

References kaldi::cblas_Xaxpy(), kaldi::cblas_Xdot(), kaldi::cblas_Xgemv(), kaldi::cblas_Xger(), kaldi::cblas_Xspmv(), kaldi::cblas_Xspr2(), VectorBase< Real >::Data(), MatrixBase< Real >::Data(), kaldi::HouseBackward(), rnnlm::i, KALDI_ASSERT, kaldi::kTrans, rnnlm::n, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), MatrixBase< Real >::SetUnit(), and MatrixBase< Real >::Stride().

Referenced by SpMatrix< Real >::Eig(), kaldi::UnitTestTridiagonalize(), and kaldi::UnitTestTridiagonalizeAndQr().

165  {
166  MatrixIndexT n = this->NumRows();
167  KALDI_ASSERT(Q == NULL || (Q->NumRows() == n &&
168  Q->NumCols() == n));
169  if (Q != NULL) Q->SetUnit();
170  Real *data = this->Data();
171  Real *qdata = (Q == NULL ? NULL : Q->Data());
172  MatrixIndexT qstride = (Q == NULL ? 0 : Q->Stride());
173  Vector<Real> tmp_v(n-1), tmp_p(n);
174  Real beta, *v = tmp_v.Data(), *p = tmp_p.Data(), *w = p, *x = p;
175  for (MatrixIndexT k = n-1; k >= 2; k--) {
176  MatrixIndexT ksize = ((k+1)*k)/2;
177  // ksize is the packed size of the lower-triangular matrix of size k,
178  // which is the size of "all rows previous to this one."
179  Real *Arow = data + ksize; // In Golub+Van Loan it was A(k+1:n, k), we
180  // have Arow = A(k, 0:k-1).
181  HouseBackward(k, Arow, v, &beta); // sets v and beta.
182  cblas_Xspmv(k, beta, data, v, 1, 0.0, p, 1); // p = beta * A(0:k-1,0:k-1) v
183  Real minus_half_beta_pv = -0.5 * beta * cblas_Xdot(k, p, 1, v, 1);
184  cblas_Xaxpy(k, minus_half_beta_pv, v, 1, w, 1); // w = p - (beta p^T v/2) v;
185  // this relies on the fact that w and p are the same pointer.
186  // We're doing A(k, k-1) = ||Arow||. It happens that this element
187  // is indexed at ksize + k - 1 in the packed lower-triangular format.
188  data[ksize + k - 1] = std::sqrt(cblas_Xdot(k, Arow, 1, Arow, 1));
189  for (MatrixIndexT i = 0; i + 1 < k; i++)
190  data[ksize + i] = 0; // This is not in Golub and Van Loan but is
191  // necessary if we're not using parts of A to store the Householder
192  // vectors.
193  // We're doing A(0:k-1,0:k-1) -= (v w' + w v')
194  cblas_Xspr2(k, -1.0, v, 1, w, 1, data);
195  if (Q != NULL) { // C.f. Golub, Q is H_1 .. H_n-2... in this
196  // case we apply them in the opposite order so it's H_n-1 .. H_1,
197  // but also Q is transposed so we really have Q = H_1 .. H_n-1.
198  // It's a double negative.
199  // Anyway, we left-multiply Q by each one. The H_n would each be
200  // diag(I + beta v v', I) but we don't ever touch the last dims.
201  // We do (in Matlab notation):
202  // Q(0:k-1,:) = (I - beta v v') * Q, i.e.:
203  // Q(:,0:i-1) += -beta v (v' Q(:,0:k-1)v .. let x = -beta Q(0:k-1,:)^T v.
204  cblas_Xgemv(kTrans, k, n, -beta, qdata, qstride, v, 1, 0.0, x, 1);
205  // now x = -beta Q(:,0:k-1) v.
206  // The next line does: Q(:,0:k-1) += v x'.
207  cblas_Xger(k, n, 1.0, v, 1, x, 1, qdata, qstride);
208  }
209  }
210 }
void HouseBackward(MatrixIndexT dim, const Real *x, Real *v, Real *beta)
Definition: qr.cc:102
void cblas_Xspr2(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, const float *Ydata, MatrixIndexT incY, float *Adata)
MatrixIndexT NumRows() const
float cblas_Xdot(const int N, const float *const X, const int incX, const float *const Y, const int incY)
int32 MatrixIndexT
Definition: matrix-common.h:96
struct rnnlm::@11::@12 n
void cblas_Xaxpy(const int N, const float alpha, const float *X, const int incX, float *Y, const int incY)
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
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)
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 cblas_Xger(MatrixIndexT num_rows, MatrixIndexT num_cols, float alpha, const float *xdata, MatrixIndexT incX, const float *ydata, MatrixIndexT incY, float *Mdata, MatrixIndexT stride)

Friends And Related Function Documentation

friend class CuSpMatrix< Real >
friend

Definition at line 42 of file sp-matrix.h.

friend class std::vector< Matrix< Real > >
friend

Definition at line 45 of file sp-matrix.h.


The documentation for this class was generated from the following files: