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...
 
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:112

Member Function Documentation

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

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

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

183  {
184  int32 num_rows = this->num_rows_;
185  KALDI_ASSERT(num_rows == v.Dim() && num_rows > 0);
186  const OtherReal *src = v.Data();
187  Real *dst = this->data_;
188  if (alpha == 1.0)
189  for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
190  *dst += *src;
191  else
192  for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
193  *dst += alpha * *src;
194 }
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 1110 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(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), OnlinePreconditionerSimple::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::UnitTestTopEigs(), and kaldi::UnitTestTp2().

1111  {
1112  KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows())
1113  || (transM == kTrans && this->NumRows() == M.NumCols()));
1114 
1115  // Cblas has no function *sprk (i.e. symmetric packed rank-k update), so we
1116  // use as temporary storage a regular matrix of which we only access its lower
1117  // triangle
1118 
1119  MatrixIndexT this_dim = this->NumRows(),
1120  m_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows());
1121 
1122  if (this_dim == 0) return;
1123  if (alpha == 0.0) {
1124  if (beta != 1.0) this->Scale(beta);
1125  return;
1126  }
1127 
1128  Matrix<Real> temp_mat(*this); // wastefully copies upper triangle too, but this
1129  // doesn't dominate O(N) time.
1130 
1131  // This function call is hard-coded to update the lower triangle.
1132  cblas_Xsyrk(transM, this_dim, m_other_dim, alpha, M.Data(),
1133  M.Stride(), beta, temp_mat.Data(), temp_mat.Stride());
1134 
1135  this->CopyFromMat(temp_mat, kTakeLower);
1136 }
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:112
#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 982 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().

984  {
985  if (transM == kNoTrans) {
986  KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
987  } else {
988  KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
989  }
990  Vector<Real> tmp_vec(A.NumRows());
991  Real *tmp_vec_data = tmp_vec.Data();
992  SpMatrix<Real> tmp_A;
993  const Real *p_A_data = A.Data();
994  Real *p_row_data = this->Data();
995  MatrixIndexT M_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows()),
996  M_same_dim = (transM == kNoTrans ? M.NumRows() : M.NumCols()),
997  M_stride = M.Stride(), dim = this->NumRows();
998  KALDI_ASSERT(M_same_dim == dim);
999 
1000  const Real *M_data = M.Data();
1001 
1002  if (this->Data() <= A.Data() + A.SizeInBytes() &&
1003  this->Data() + this->SizeInBytes() >= A.Data()) {
1004  // Matrices A and *this overlap. Make copy of A
1005  tmp_A.Resize(A.NumRows());
1006  tmp_A.CopyFromSp(A);
1007  p_A_data = tmp_A.Data();
1008  }
1009 
1010  if (transM == kNoTrans) {
1011  for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
1012  cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.RowData(r), 1, 0.0, tmp_vec_data, 1);
1013  cblas_Xgemv(transM, r+1, M_other_dim, alpha, M_data, M_stride,
1014  tmp_vec_data, 1, beta, p_row_data, 1);
1015  }
1016  } else {
1017  for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
1018  cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.Data() + r, M.Stride(), 0.0, tmp_vec_data, 1);
1019  cblas_Xgemv(transM, M_other_dim, r+1, alpha, M_data, M_stride,
1020  tmp_vec_data, 1, beta, p_row_data, 1);
1021  }
1022  }
1023 }
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 1081 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().

1085  {
1086  this->Scale(beta);
1087  KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows() &&
1088  M.NumCols() == v.Dim()) ||
1089  (transM == kTrans && this->NumRows() == M.NumCols() &&
1090  M.NumRows() == v.Dim()));
1091 
1092  if (transM == kNoTrans) {
1093  const Real *Mdata = M.Data(), *vdata = v.Data();
1094  Real *data = this->data_;
1095  MatrixIndexT dim = this->NumRows(), mcols = M.NumCols(),
1096  mstride = M.Stride();
1097  for (MatrixIndexT col = 0; col < mcols; col++, vdata++, Mdata += 1)
1098  cblas_Xspr(dim, *vdata*alpha, Mdata, mstride, data);
1099  } else {
1100  const Real *Mdata = M.Data(), *vdata = v.Data();
1101  Real *data = this->data_;
1102  MatrixIndexT dim = this->NumRows(), mrows = M.NumRows(),
1103  mstride = M.Stride();
1104  for (MatrixIndexT row = 0; row < mrows; row++, vdata++, Mdata += mstride)
1105  cblas_Xspr(dim, *vdata*alpha, Mdata, 1, data);
1106  }
1107 }
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 1026 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().

1029  {
1030  KALDI_ASSERT((transM == kNoTrans && M.NumCols() == A.NumRows()) ||
1031  (transM == kTrans && M.NumRows() == A.NumRows()));
1032  if (transM == kNoTrans) {
1033  KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
1034  } else {
1035  KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
1036  }
1037  MatrixIndexT Adim = A.NumRows(), dim = this->num_rows_;
1038 
1039  Matrix<Real> temp_A(A); // represent A as full matrix.
1040  Matrix<Real> temp_MA(dim, Adim);
1041  temp_MA.AddSmatMat(1.0, M, transM, temp_A, kNoTrans, 0.0);
1042 
1043  // Next-- we want to do *this = alpha * temp_MA * M^T + beta * *this.
1044  // To make it sparse vector multiplies, since M is sparse, we'd like
1045  // to do: for each column c, (*this column c) += temp_MA * (M^T's column c.)
1046  // [ignoring the alpha and beta here.]
1047  // It's not convenient to process columns in the symmetric
1048  // packed format because they don't have a constant stride. However,
1049  // we can use the fact that temp_MA * M is symmetric, to just assign
1050  // each row of *this instead of each column.
1051  // So the final iteration is:
1052  // for i = 0... dim-1,
1053  // [the i'th row of *this] = beta * [the i'th row of *this] + alpha *
1054  // temp_MA * [the i'th column of M].
1055  // Of course, we only process the first 0 ... i elements of this row,
1056  // as that's all that are kept in the symmetric packed format.
1057 
1058  Matrix<Real> temp_this(*this);
1059  Real *data = this->data_;
1060  const Real *Mdata = M.Data(), *MAdata = temp_MA.Data();
1061  MatrixIndexT temp_MA_stride = temp_MA.Stride(), Mstride = M.Stride();
1062 
1063  if (transM == kNoTrans) {
1064  // The column of M^T corresponds to the rows of the supplied matrix.
1065  for (MatrixIndexT i = 0; i < dim; i++, data += i) {
1066  MatrixIndexT num_rows = i + 1, num_cols = Adim;
1067  Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
1068  temp_MA_stride, Mdata + (i * Mstride), 1, beta, data, 1);
1069  }
1070  } else {
1071  // The column of M^T corresponds to the columns of the supplied matrix.
1072  for (MatrixIndexT i = 0; i < dim; i++, data += i) {
1073  MatrixIndexT num_rows = i + 1, num_cols = Adim;
1074  Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
1075  temp_MA_stride, Mdata + i, Mstride, beta, data, 1);
1076  }
1077  }
1078 }
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 211 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().

211  {
212  this->AddPacked(alpha, Ma);
213  }
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 1156 of file sp-matrix.cc.

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

1157  {
1158  Matrix<Real> Tmat(T);
1159  AddMat2(alpha, Tmat, transM, beta);
1160 }
void AddMat2(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const Real beta)
rank-N update: if (transM == kNoTrans) (*this) = beta*(*this) + alpha * M * M^T, or (if transM == kTr...
Definition: sp-matrix.cc:1110
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 1139 of file sp-matrix.cc.

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

1141  {
1142  Matrix<Real> Tmat(T);
1143  AddMat2Sp(alpha, Tmat, transM, A, beta);
1144 }
void AddMat2Sp(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const SpMatrix< Real > &A, const Real beta=0.0)
Extension of rank-N update: this <– beta*this + alpha * M * A * M^T.
Definition: sp-matrix.cc:982
void AddVec2 ( const double  alpha,
const VectorBase< double > &  v 
)
void AddVec2 ( const Real  alpha,
const VectorBase< OtherReal > &  v 
)

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

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

946  {
947  KALDI_ASSERT(v.Dim() == this->NumRows());
948  Real *data = this->data_;
949  const OtherReal *v_data = v.Data();
950  MatrixIndexT nr = this->num_rows_;
951  for (MatrixIndexT i = 0; i < nr; i++)
952  for (MatrixIndexT j = 0; j <= i; j++, data++)
953  *data += alpha * v_data[i] * v_data[j];
954 }
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 float  alpha,
const VectorBase< float > &  v 
)

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

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

915  {
916  KALDI_ASSERT(v.Dim() == this->NumRows());
917  cblas_Xspr(v.Dim(), alpha, v.Data(), 1,
918  this->data_);
919 }
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 938 of file sp-matrix.cc.

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

938  {
939  KALDI_ASSERT(v.Dim() == num_rows_);
940  cblas_Xspr(v.Dim(), alpha, v.Data(), 1, data_);
941 }
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 922 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().

923  {
924  KALDI_ASSERT(v.Dim() == this->NumRows() && S.NumRows() == this->NumRows());
925  const Real *Sdata = S.Data();
926  const Real *vdata = v.Data();
927  Real *data = this->data_;
928  MatrixIndexT dim = this->num_rows_;
929  for (MatrixIndexT r = 0; r < dim; r++)
930  for (MatrixIndexT c = 0; c <= r; c++, Sdata++, data++)
931  *data = beta * *data + alpha * vdata[r] * vdata[c] * *Sdata;
932 }
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 1147 of file sp-matrix.cc.

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

Referenced by kaldi::UnitTestSpAddVecVec().

1148  {
1149  int32 dim = this->NumRows();
1150  KALDI_ASSERT(dim == v.Dim() && dim == w.Dim() && dim > 0);
1151  cblas_Xspr2(dim, alpha, v.Data(), 1, w.Data(), 1, this->data_);
1152 }
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 536 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().

537  {
538  MatrixIndexT dim = this->NumRows();
539  int nfloored = 0;
540  KALDI_ASSERT(C.NumRows() == dim);
541  KALDI_ASSERT(alpha > 0);
542  TpMatrix<Real> L(dim);
543  L.Cholesky(C);
544  L.Scale(std::sqrt(alpha)); // equivalent to scaling C by alpha.
545  TpMatrix<Real> LInv(L);
546  LInv.Invert();
547 
548  SpMatrix<Real> D(dim);
549  { // D = L^{-1} * (*this) * L^{-T}
550  Matrix<Real> LInvFull(LInv);
551  D.AddMat2Sp(1.0, LInvFull, kNoTrans, (*this), 0.0);
552  }
553 
554  Vector<Real> l(dim);
555  Matrix<Real> U(dim, dim);
556 
557  D.Eig(&l, &U);
558 
559  if (verbose) {
560  KALDI_LOG << "ApplyFloor: flooring following diagonal to 1: " << l;
561  }
562  for (MatrixIndexT i = 0; i < l.Dim(); i++) {
563  if (l(i) < 1.0) {
564  nfloored++;
565  l(i) = 1.0;
566  }
567  }
568  l.ApplyPow(0.5);
569  U.MulColsVec(l);
570  D.AddMat2(1.0, U, kNoTrans, 0.0);
571  { // D' := U * diag(l') * U^T ... l'=floor(l, 1)
572  Matrix<Real> LFull(L);
573  (*this).AddMat2Sp(1.0, LFull, kNoTrans, D, 0.0); // A := L * D' * L^T
574  }
575  return nfloored;
576 }
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 589 of file sp-matrix.cc.

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

589  {
590  MatrixIndexT Dim = this->NumRows();
591  int nfloored = 0;
592  Vector<Real> s(Dim);
593  Matrix<Real> P(Dim, Dim);
594  (*this).Eig(&s, &P);
595  for (MatrixIndexT i = 0; i < Dim; i++) {
596  if (s(i) < floor) {
597  nfloored++;
598  s(i) = floor;
599  }
600  }
601  (*this).AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
602  return nfloored;
603 }
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 91 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().

91  {
92  if (power == 1) return; // can do nothing.
93  MatrixIndexT D = this->NumRows();
94  KALDI_ASSERT(D > 0);
95  Matrix<Real> U(D, D);
96  Vector<Real> l(D);
97  (*this).SymPosSemiDefEig(&l, &U);
98 
99  Vector<Real> l_copy(l);
100  try {
101  l.ApplyPow(power * 0.5);
102  }
103  catch(...) {
104  KALDI_ERR << "Error taking power " << (power * 0.5) << " of vector "
105  << l_copy;
106  }
107  U.MulColsVec(l);
108  (*this).AddMat2(1.0, U, kNoTrans, 0.0);
109 }
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 525 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().

525  {
526  if (this->NumRows() != other.NumRows())
527  KALDI_ERR << "SpMatrix::AproxEqual, size mismatch, "
528  << this->NumRows() << " vs. " << other.NumRows();
529  SpMatrix<Real> tmp(*this);
530  tmp.AddSp(-1.0, other);
531  return (tmp.FrobeniusNorm() <= tol * std::max(this->FrobeniusNorm(), other.FrobeniusNorm()));
532 }
Real FrobeniusNorm() const
sqrt of sum of square elements.
Definition: sp-matrix.cc:513
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 112 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 BasisFmllrEstimate::ComputeAmDiagPrecond(), SpMatrix< BaseFloat >::SpMatrix(), kaldi::UnitTestCopySp(), kaldi::UnitTestDeterminant(), kaldi::UnitTestDeterminantSign(), and kaldi::UnitTestPower().

113  {
114  KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
115  MatrixIndexT D = this->NumRows();
116 
117  switch (copy_type) {
118  case kTakeMeanAndCheck:
119  {
120  Real good_sum = 0.0, bad_sum = 0.0;
121  for (MatrixIndexT i = 0; i < D; i++) {
122  for (MatrixIndexT j = 0; j < i; j++) {
123  Real a = M(i, j), b = M(j, i), avg = 0.5*(a+b), diff = 0.5*(a-b);
124  (*this)(i, j) = avg;
125  good_sum += std::abs(avg);
126  bad_sum += std::abs(diff);
127  }
128  good_sum += std::abs(M(i, i));
129  (*this)(i, i) = M(i, i);
130  }
131  if (bad_sum > 0.01 * good_sum) {
132  KALDI_ERR << "SpMatrix::Copy(), source matrix is not symmetric: "
133  << bad_sum << ">" << good_sum;
134  }
135  break;
136  }
137  case kTakeMean:
138  {
139  for (MatrixIndexT i = 0; i < D; i++) {
140  for (MatrixIndexT j = 0; j < i; j++) {
141  (*this)(i, j) = 0.5*(M(i, j) + M(j, i));
142  }
143  (*this)(i, i) = M(i, i);
144  }
145  break;
146  }
147  case kTakeLower:
148  { // making this one a bit more efficient.
149  const Real *src = M.Data();
150  Real *dest = this->data_;
151  MatrixIndexT stride = M.Stride();
152  for (MatrixIndexT i = 0; i < D; i++) {
153  for (MatrixIndexT j = 0; j <= i; j++)
154  dest[j] = src[j];
155  dest += i + 1;
156  src += stride;
157  }
158  }
159  break;
160  case kTakeUpper:
161  for (MatrixIndexT i = 0; i < D; i++)
162  for (MatrixIndexT j = 0; j <= i; j++)
163  (*this)(i, j) = M(j, i);
164  break;
165  default:
166  KALDI_ASSERT("Invalid argument to SpMatrix::CopyFromMat");
167  }
168 }
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
Real FrobeniusNorm ( ) const

sqrt of sum of square elements.

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

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

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

513  {
514  Real sum = 0.0;
515  MatrixIndexT R = this->NumRows();
516  for (MatrixIndexT i = 0; i < R; i++) {
517  for (MatrixIndexT j = 0; j < i; j++)
518  sum += (*this)(i, j) * (*this)(i, j) * 2;
519  sum += (*this)(i, i) * (*this)(i, i);
520  }
521  return std::sqrt(sum);
522 }
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 219 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().

219  {
220  // these are CLAPACK types
221  KaldiBlasInt result;
222  KaldiBlasInt rows = static_cast<int>(this->num_rows_);
223  KaldiBlasInt* p_ipiv = new KaldiBlasInt[rows];
224  Real *p_work; // workspace for the lapack function
225  void *temp;
226  if ((p_work = static_cast<Real*>(
227  KALDI_MEMALIGN(16, sizeof(Real) * rows, &temp))) == NULL) {
228  delete[] p_ipiv;
229  throw std::bad_alloc();
230  }
231 #ifdef HAVE_OPENBLAS
232  memset(p_work, 0, sizeof(Real) * rows); // gets rid of a probably
233  // spurious Valgrind warning about jumps depending upon uninitialized values.
234 #endif
235 
236 
237  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
238  // of the data. We have a row-wise storage, therefore, we need to "invert"
239  clapack_Xsptrf(&rows, this->data_, p_ipiv, &result);
240 
241 
242  KALDI_ASSERT(result >= 0 && "Call to CLAPACK ssptrf_ called with wrong arguments");
243 
244  if (result > 0) { // Singular...
245  if (det_sign) *det_sign = 0;
246  if (logdet) *logdet = -std::numeric_limits<Real>::infinity();
247  if (need_inverse) KALDI_ERR << "CLAPACK stptrf_ : factorization failed";
248  } else { // Not singular.. compute log-determinant if needed.
249  if (logdet != NULL || det_sign != NULL) {
250  Real prod = 1.0, log_prod = 0.0;
251  int sign = 1;
252  for (int i = 0; i < (int)this->num_rows_; i++) {
253  if (p_ipiv[i] > 0) { // not a 2x2 block...
254  // if (p_ipiv[i] != i+1) sign *= -1; // row swap.
255  Real diag = (*this)(i, i);
256  prod *= diag;
257  } else { // negative: 2x2 block. [we are in first of the two].
258  i++; // skip over the first of the pair.
259  // each 2x2 block...
260  Real diag1 = (*this)(i, i), diag2 = (*this)(i-1, i-1),
261  offdiag = (*this)(i, i-1);
262  Real thisdet = diag1*diag2 - offdiag*offdiag;
263  // thisdet == determinant of 2x2 block.
264  // The following line is more complex than it looks: there are 2 offsets of
265  // 1 that cancel.
266  prod *= thisdet;
267  }
268  if (i == (int)(this->num_rows_-1) || fabs(prod) < 1.0e-10 || fabs(prod) > 1.0e+10) {
269  if (prod < 0) { prod = -prod; sign *= -1; }
270  log_prod += kaldi::Log(std::abs(prod));
271  prod = 1.0;
272  }
273  }
274  if (logdet != NULL) *logdet = log_prod;
275  if (det_sign != NULL) *det_sign = sign;
276  }
277  }
278  if (!need_inverse) {
279  delete [] p_ipiv;
280  KALDI_MEMALIGN_FREE(p_work);
281  return; // Don't need what is computed next.
282  }
283  // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
284  // of the data. We have a row-wise storage, therefore, we need to "invert"
285  clapack_Xsptri(&rows, this->data_, p_ipiv, p_work, &result);
286 
287  KALDI_ASSERT(result >=0 &&
288  "Call to CLAPACK ssptri_ called with wrong arguments");
289 
290  if (result != 0) {
291  KALDI_ERR << "CLAPACK ssptrf_ : Matrix is singular";
292  }
293 
294  delete [] p_ipiv;
295  KALDI_MEMALIGN_FREE(p_work);
296 }
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 312 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().

313  {
314  SpMatrix<double> dmat(*this);
315  double logdet_tmp, det_sign_tmp;
316  dmat.Invert(logdet ? &logdet_tmp : NULL,
317  det_sign ? &det_sign_tmp : NULL,
318  inverse_needed);
319  if (logdet) *logdet = logdet_tmp;
320  if (det_sign) *det_sign = det_sign_tmp;
321  (*this).CopyFromSp(dmat);
322 }
bool IsDiagonal ( Real  cutoff = 1.0e-05) const

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

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

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

465  {
466  MatrixIndexT R = this->NumRows();
467  Real bad_sum = 0.0, good_sum = 0.0;
468  for (MatrixIndexT i = 0; i < R; i++) {
469  for (MatrixIndexT j = 0; j <= i; j++) {
470  if (i == j)
471  good_sum += std::abs((*this)(i, j));
472  else
473  bad_sum += std::abs((*this)(i, j));
474  }
475  }
476  return (!(bad_sum > good_sum * cutoff));
477 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
bool IsPosDef ( ) const

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

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

Referenced by MleAmSgmm2Updater::RenormalizeV().

75  {
76  MatrixIndexT D = (*this).NumRows();
77  KALDI_ASSERT(D > 0);
78  try {
79  TpMatrix<Real> C(D);
80  C.Cholesky(*this);
81  for (MatrixIndexT r = 0; r < D; r++)
82  if (C(r, r) == 0.0) return false;
83  return true;
84  }
85  catch(...) { // not positive semidefinite.
86  return false;
87  }
88 }
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 491 of file sp-matrix.cc.

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

Referenced by kaldi::UnitTestTridiag().

491  {
492  MatrixIndexT R = this->NumRows();
493  Real max_abs_2diag = 0.0, max_abs_offdiag = 0.0;
494  for (MatrixIndexT i = 0; i < R; i++)
495  for (MatrixIndexT j = 0; j <= i; j++) {
496  if (j+1 < i)
497  max_abs_offdiag = std::max(max_abs_offdiag,
498  std::abs((*this)(i, j)));
499  else
500  max_abs_2diag = std::max(max_abs_2diag,
501  std::abs((*this)(i, j)));
502  }
503  return (max_abs_offdiag <= cutoff * max_abs_2diag);
504 }
MatrixIndexT NumRows() const
int32 MatrixIndexT
Definition: matrix-common.h:96
bool IsZero ( Real  cutoff = 1.0e-05) const

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

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

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

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

606  { // e.g. maxCond = 1.0e+05.
607  MatrixIndexT Dim = this->NumRows();
608  Vector<Real> s(Dim);
609  Matrix<Real> P(Dim, Dim);
610  (*this).SymPosSemiDefEig(&s, &P);
611  KALDI_ASSERT(maxCond > 1);
612  Real floor = s.Max() / maxCond;
613  if (floor < 0) floor = 0;
614  if (floor < 1.0e-40) {
615  KALDI_WARN << "LimitCond: limiting " << floor << " to 1.0e-40";
616  floor = 1.0e-40;
617  }
618  MatrixIndexT nfloored = 0;
619  for (MatrixIndexT i = 0; i < Dim; i++) {
620  if (s(i) <= floor) nfloored++;
621  if (invert)
622  s(i) = 1.0 / std::sqrt(std::max(s(i), floor));
623  else
624  s(i) = std::sqrt(std::max(s(i), floor));
625  }
626  P.MulColsVec(s);
627  (*this).AddMat2(1.0, P, kNoTrans, 0.0); // (*this) = P*P^T. ... (*this) = P * floor(s) * P^T ... if P was original P.
628  return nfloored;
629 }
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 333 of file sp-matrix.h.

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

333  {
334  SpMatrix<double> dmat(*this);
335  MatrixIndexT ans = dmat.LimitCond(maxCond, invert);
336  (*this).CopyFromSp(dmat);
337  return ans;
338  }
int32 MatrixIndexT
Definition: matrix-common.h:96
Real LogDet ( Real *  det_sign = NULL) const

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

References SpMatrix< Real >::Invert().

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

579  {
580  Real log_det;
581  SpMatrix<Real> tmp(*this);
582  // false== output not needed (saves some computation).
583  tmp.Invert(&log_det, det_sign, false);
584  return log_det;
585 }
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 203 of file sp-matrix.h.

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

203  {
204  Vector<Real> s((*this).NumRows());
205  Matrix<Real> P((*this).NumRows(), (*this).NumCols());
206  SymPosSemiDefEig(&s, &P);
207  KALDI_LOG << "PrintEigs: " << name << ": " << s;
208  }
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:491
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: