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

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

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

1111  {
1112  this->Scale(beta);
1113  KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows() &&
1114  M.NumCols() == v.Dim()) ||
1115  (transM == kTrans && this->NumRows() == M.NumCols() &&
1116  M.NumRows() == v.Dim()));
1117 
1118  if (transM == kNoTrans) {
1119  const Real *Mdata = M.Data(), *vdata = v.Data();
1120  Real *data = this->data_;
1121  MatrixIndexT dim = this->NumRows(), mcols = M.NumCols(),
1122  mstride = M.Stride();
1123  for (MatrixIndexT col = 0; col < mcols; col++, vdata++, Mdata += 1)
1124  cblas_Xspr(dim, *vdata*alpha, Mdata, mstride, data);
1125  } else {
1126  const Real *Mdata = M.Data(), *vdata = v.Data();
1127  Real *data = this->data_;
1128  MatrixIndexT dim = this->NumRows(), mrows = M.NumRows(),
1129  mstride = M.Stride();
1130  for (MatrixIndexT row = 0; row < mrows; row++, vdata++, Mdata += mstride)
1131  cblas_Xspr(dim, *vdata*alpha, Mdata, 1, data);
1132  }
1133 }
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 1052 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().

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

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

1183  {
1184  Matrix<Real> Tmat(T);
1185  AddMat2(alpha, Tmat, transM, beta);
1186 }
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:1136
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 1165 of file sp-matrix.cc.

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

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

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

Definition at line 972 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().

972  {
973  KALDI_ASSERT(v.Dim() == this->NumRows());
974  Real *data = this->data_;
975  const OtherReal *v_data = v.Data();
976  MatrixIndexT nr = this->num_rows_;
977  for (MatrixIndexT i = 0; i < nr; i++)
978  for (MatrixIndexT j = 0; j <= i; j++, data++)
979  *data += alpha * v_data[i] * v_data[j];
980 }
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 941 of file sp-matrix.cc.

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

941  {
942  KALDI_ASSERT(v.Dim() == this->NumRows());
943  cblas_Xspr(v.Dim(), alpha, v.Data(), 1,
944  this->data_);
945 }
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 964 of file sp-matrix.cc.

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

964  {
965  KALDI_ASSERT(v.Dim() == num_rows_);
966  cblas_Xspr(v.Dim(), alpha, v.Data(), 1, data_);
967 }
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 948 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().

949  {
950  KALDI_ASSERT(v.Dim() == this->NumRows() && S.NumRows() == this->NumRows());
951  const Real *Sdata = S.Data();
952  const Real *vdata = v.Data();
953  Real *data = this->data_;
954  MatrixIndexT dim = this->num_rows_;
955  for (MatrixIndexT r = 0; r < dim; r++)
956  for (MatrixIndexT c = 0; c <= r; c++, Sdata++, data++)
957  *data = beta * *data + alpha * vdata[r] * vdata[c] * *Sdata;
958 }
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 1173 of file sp-matrix.cc.

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

Referenced by kaldi::UnitTestSpAddVecVec().

1174  {
1175  int32 dim = this->NumRows();
1176  KALDI_ASSERT(dim == v.Dim() && dim == w.Dim() && dim > 0);
1177  cblas_Xspr2(dim, alpha, v.Data(), 1, w.Data(), 1, this->data_);
1178 }
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 562 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().

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

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

615  {
616  MatrixIndexT Dim = this->NumRows();
617  int nfloored = 0;
618  Vector<Real> s(Dim);
619  Matrix<Real> P(Dim, Dim);
620  (*this).Eig(&s, &P);
621  for (MatrixIndexT i = 0; i < Dim; i++) {
622  if (s(i) < floor) {
623  nfloored++;
624  s(i) = floor;
625  }
626  }
627  (*this).AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
628  return nfloored;
629 }
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 551 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().

551  {
552  if (this->NumRows() != other.NumRows())
553  KALDI_ERR << "SpMatrix::AproxEqual, size mismatch, "
554  << this->NumRows() << " vs. " << other.NumRows();
555  SpMatrix<Real> tmp(*this);
556  tmp.AddSp(-1.0, other);
557  return (tmp.FrobeniusNorm() <= tol * std::max(this->FrobeniusNorm(), other.FrobeniusNorm()));
558 }
Real FrobeniusNorm() const
sqrt of sum of square elements.
Definition: sp-matrix.cc:539
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(), OnlinePreconditioner::PreconditionDirectionsInternal(), OnlineNaturalGradient::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 539 of file sp-matrix.cc.

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

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

539  {
540  Real sum = 0.0;
541  MatrixIndexT R = this->NumRows();
542  for (MatrixIndexT i = 0; i < R; i++) {
543  for (MatrixIndexT j = 0; j < i; j++)
544  sum += (*this)(i, j) * (*this)(i, j) * 2;
545  sum += (*this)(i, i) * (*this)(i, i);
546  }
547  return std::sqrt(sum);
548 }
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  delete[] p_ipiv;
255  throw std::bad_alloc();
256  }
257 #ifdef HAVE_OPENBLAS
258  memset(p_work, 0, sizeof(Real) * rows); // gets rid of a probably
259  // spurious Valgrind warning about jumps depending upon uninitialized values.
260 #endif
261 
262 
263  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
264  // of the data. We have a row-wise storage, therefore, we need to "invert"
265  clapack_Xsptrf(&rows, this->data_, p_ipiv, &result);
266 
267 
268  KALDI_ASSERT(result >= 0 && "Call to CLAPACK ssptrf_ called with wrong arguments");
269 
270  if (result > 0) { // Singular...
271  if (det_sign) *det_sign = 0;
272  if (logdet) *logdet = -std::numeric_limits<Real>::infinity();
273  if (need_inverse) KALDI_ERR << "CLAPACK stptrf_ : factorization failed";
274  } else { // Not singular.. compute log-determinant if needed.
275  if (logdet != NULL || det_sign != NULL) {
276  Real prod = 1.0, log_prod = 0.0;
277  int sign = 1;
278  for (int i = 0; i < (int)this->num_rows_; i++) {
279  if (p_ipiv[i] > 0) { // not a 2x2 block...
280  // if (p_ipiv[i] != i+1) sign *= -1; // row swap.
281  Real diag = (*this)(i, i);
282  prod *= diag;
283  } else { // negative: 2x2 block. [we are in first of the two].
284  i++; // skip over the first of the pair.
285  // each 2x2 block...
286  Real diag1 = (*this)(i, i), diag2 = (*this)(i-1, i-1),
287  offdiag = (*this)(i, i-1);
288  Real thisdet = diag1*diag2 - offdiag*offdiag;
289  // thisdet == determinant of 2x2 block.
290  // The following line is more complex than it looks: there are 2 offsets of
291  // 1 that cancel.
292  prod *= thisdet;
293  }
294  if (i == (int)(this->num_rows_-1) || fabs(prod) < 1.0e-10 || fabs(prod) > 1.0e+10) {
295  if (prod < 0) { prod = -prod; sign *= -1; }
296  log_prod += kaldi::Log(std::abs(prod));
297  prod = 1.0;
298  }
299  }
300  if (logdet != NULL) *logdet = log_prod;
301  if (det_sign != NULL) *det_sign = sign;
302  }
303  }
304  if (!need_inverse) {
305  delete [] p_ipiv;
306  KALDI_MEMALIGN_FREE(p_work);
307  return; // Don't need what is computed next.
308  }
309  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
310  // of the data. We have a row-wise storage, therefore, we need to "invert"
311  clapack_Xsptri(&rows, this->data_, p_ipiv, p_work, &result);
312 
313  KALDI_ASSERT(result >=0 &&
314  "Call to CLAPACK ssptri_ called with wrong arguments");
315 
316  if (result != 0) {
317  KALDI_ERR << "CLAPACK ssptrf_ : Matrix is singular";
318  }
319 
320  delete [] p_ipiv;
321  KALDI_MEMALIGN_FREE(p_work);
322 }
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 338 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().

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

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

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

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

491  {
492  MatrixIndexT R = this->NumRows();
493  Real bad_sum = 0.0, good_sum = 0.0;
494  for (MatrixIndexT i = 0; i < R; i++) {
495  for (MatrixIndexT j = 0; j <= i; j++) {
496  if (i == j)
497  good_sum += std::abs((*this)(i, j));
498  else
499  bad_sum += std::abs((*this)(i, j));
500  }
501  }
502  return (!(bad_sum > good_sum * cutoff));
503 }
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 517 of file sp-matrix.cc.

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

Referenced by kaldi::UnitTestTridiag().

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

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

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

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

Definition at line 632 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().

632  { // e.g. maxCond = 1.0e+05.
633  MatrixIndexT Dim = this->NumRows();
634  Vector<Real> s(Dim);
635  Matrix<Real> P(Dim, Dim);
636  (*this).SymPosSemiDefEig(&s, &P);
637  KALDI_ASSERT(maxCond > 1);
638  Real floor = s.Max() / maxCond;
639  if (floor < 0) floor = 0;
640  if (floor < 1.0e-40) {
641  KALDI_WARN << "LimitCond: limiting " << floor << " to 1.0e-40";
642  floor = 1.0e-40;
643  }
644  MatrixIndexT nfloored = 0;
645  for (MatrixIndexT i = 0; i < Dim; i++) {
646  if (s(i) <= floor) nfloored++;
647  if (invert)
648  s(i) = 1.0 / std::sqrt(std::max(s(i), floor));
649  else
650  s(i) = std::sqrt(std::max(s(i), floor));
651  }
652  P.MulColsVec(s);
653  (*this).AddMat2(1.0, P, kNoTrans, 0.0); // (*this) = P*P^T. ... (*this) = P * floor(s) * P^T ... if P was original P.
654  return nfloored;
655 }
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:1107
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 605 of file sp-matrix.cc.

References SpMatrix< Real >::Invert().

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

605  {
606  Real log_det;
607  SpMatrix<Real> tmp(*this);
608  // false== output not needed (saves some computation).
609  tmp.Invert(&log_det, det_sign, false);
610  return log_det;
611 }
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:517
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: