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 <– betathis + 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 <– betathis + 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 <– betathis + 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 <– betathis + 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 62 of file matrix-common.h.

Constructor & Destructor Documentation

◆ SpMatrix() [1/6]

SpMatrix ( )

inline

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

Referenced by SpMatrix< float >::SpMatrix().

47 : PackedMatrix<Real>() {}

◆ SpMatrix() [2/6]

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.

                                                    {
    Resize(cu.NumRows());
    cu.CopyToSp(this);
 }

◆ SpMatrix() [3/6]

SpMatrix	(	MatrixIndexT	r,
		MatrixResizeType	resize_type = `kSetZero`
	)

inlineexplicit

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

55 : PackedMatrix<Real>(r, resize_type) {}

◆ SpMatrix() [4/6]

SpMatrix ( const SpMatrix< Real > & orig )

inline

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

58 : PackedMatrix<Real>(orig) {}

◆ SpMatrix() [5/6]

SpMatrix ( const SpMatrix< OtherReal > & orig )

inlineexplicit

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

62 : PackedMatrix<Real>(orig) {}

◆ SpMatrix() [6/6]

SpMatrix	(	const MatrixBase< Real > &	orig,
		SpCopyType	copy_type = `kTakeMean`
	)

inlineexplicit

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

       : PackedMatrix<Real>(orig.NumRows(), kUndefined) {
     CopyFromMat(orig, copy_type);
   }

Member Function Documentation

◆ AddDiagVec()

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.

Referenced by SpMatrix< float >::AddDiagVec(), SpMatrix< float >::AddSp(), Fmpe::ComputeC(), main(), and kaldi::UnitTestSpAddDiagVec().

                                                                                  {
   int32 num_rows = this->num_rows_;
   KALDI_ASSERT(num_rows == v.Dim() && num_rows > 0);
   const OtherReal *src = v.Data();
   Real *dst = this->data_;
   if (alpha == 1.0)
     for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
       *dst += *src;
   else
     for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
       *dst += alpha * *src;
 }

◆ AddMat2()

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.

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

                                                                            {
   KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows())
                || (transM == kTrans && this->NumRows() == M.NumCols()));
 
   // Cblas has no function *sprk (i.e. symmetric packed rank-k update), so we
   // use as temporary storage a regular matrix of which we only access its lower
   // triangle
 
   MatrixIndexT this_dim = this->NumRows(),
       m_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows());
 
   if (this_dim == 0) return;
   if (alpha == 0.0) {
     if (beta != 1.0) this->Scale(beta);
     return;
   }
 
   Matrix<Real> temp_mat(*this); // wastefully copies upper triangle too, but this
   // doesn't dominate O(N) time.
 
   // This function call is hard-coded to update the lower triangle.
   cblas_Xsyrk(transM, this_dim, m_other_dim, alpha, M.Data(),
               M.Stride(), beta, temp_mat.Data(), temp_mat.Stride());
 
   this->CopyFromMat(temp_mat, kTakeLower);
 }

◆ AddMat2Sp()

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.

Referenced by SpMatrix< float >::AddSp(), kaldi::ApplyFeatureTransformToStats(), SpMatrix< float >::ApplyFloor(), Plda::ApplyTransform(), 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().

                                                                           {
   if (transM == kNoTrans) {
     KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
   } else {
     KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
   }
   Vector<Real> tmp_vec(A.NumRows());
   Real *tmp_vec_data = tmp_vec.Data();
   SpMatrix<Real> tmp_A;
   const Real *p_A_data = A.Data();
   Real *p_row_data = this->Data();
   MatrixIndexT M_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows()),
       M_same_dim = (transM == kNoTrans ? M.NumRows() : M.NumCols()),
       M_stride = M.Stride(), dim = this->NumRows();
   KALDI_ASSERT(M_same_dim == dim);
 
   const Real *M_data = M.Data();
 
   if (this->Data() <= A.Data() + A.SizeInBytes() &&
       this->Data() + this->SizeInBytes() >= A.Data()) {
     // Matrices A and *this overlap. Make copy of A
     tmp_A.Resize(A.NumRows());
     tmp_A.CopyFromSp(A);
     p_A_data = tmp_A.Data();
   }
 
   if (transM == kNoTrans) {
     for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
       cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.RowData(r), 1, 0.0, tmp_vec_data, 1);
       cblas_Xgemv(transM, r+1, M_other_dim, alpha, M_data, M_stride,
                   tmp_vec_data, 1, beta, p_row_data, 1);
     }
   } else {
     for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
       cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.Data() + r, M.Stride(), 0.0, tmp_vec_data, 1);
       cblas_Xgemv(transM, M_other_dim, r+1, alpha, M_data, M_stride,
                   tmp_vec_data, 1, beta, p_row_data, 1);
     }
   }
 }

◆ AddMat2Vec()

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.

Referenced by SpMatrix< float >::AddSp(), IvectorExtractor::GetAcousticAuxfWeight(), IvectorExtractor::GetIvectorDistWeight(), IvectorExtractor::InvertWithFlooring(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), kaldi::UnitTestEigSp(), and kaldi::UnitTestMat2Vec().

                                                  {
   this->Scale(beta);
   KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows() &&
                 M.NumCols() == v.Dim()) ||
                (transM == kTrans && this->NumRows() == M.NumCols() &&
                 M.NumRows() == v.Dim()));
 
   if (transM == kNoTrans) {
     const Real *Mdata = M.Data(), *vdata = v.Data();
     Real *data = this->data_;
     MatrixIndexT dim = this->NumRows(), mcols = M.NumCols(),
         mstride = M.Stride();
     for (MatrixIndexT col = 0; col < mcols; col++, vdata++, Mdata += 1)
       cblas_Xspr(dim, *vdata*alpha, Mdata, mstride, data);
   } else {
     const Real *Mdata = M.Data(), *vdata = v.Data();
     Real *data = this->data_;
     MatrixIndexT dim = this->NumRows(), mrows = M.NumRows(),
         mstride = M.Stride();
     for (MatrixIndexT row = 0; row < mrows; row++, vdata++, Mdata += mstride)
       cblas_Xspr(dim, *vdata*alpha, Mdata, 1, data);
   }
 }

◆ AddSmat2Sp()

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.

Referenced by SpMatrix< float >::AddSp(), FmllrRawAccs::ConvertToSimpleStats(), and kaldi::UnitTestAddMat2Sp().

                      {
   KALDI_ASSERT((transM == kNoTrans && M.NumCols() == A.NumRows()) ||
                (transM == kTrans && M.NumRows() == A.NumRows()));
   if (transM == kNoTrans) {
     KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
   } else {
     KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
   }
   MatrixIndexT Adim = A.NumRows(), dim = this->num_rows_;
 
   Matrix<Real> temp_A(A); // represent A as full matrix.
   Matrix<Real> temp_MA(dim, Adim);
   temp_MA.AddSmatMat(1.0, M, transM, temp_A, kNoTrans, 0.0);
 
   // Next-- we want to do *this = alpha * temp_MA * M^T + beta * *this.
   // To make it sparse vector multiplies, since M is sparse, we'd like
   // to do: for each column c, (*this column c) += temp_MA * (M^T's column c.)
   // [ignoring the alpha and beta here.]
   // It's not convenient to process columns in the symmetric
   // packed format because they don't have a constant stride.  However,
   // we can use the fact that temp_MA * M is symmetric, to just assign
   // each row of *this instead of each column.
   // So the final iteration is:
   // for i = 0... dim-1,
   //   [the i'th row of *this] = beta * [the i'th row of *this] + alpha *
   //                               temp_MA * [the i'th column of M].
   // Of course, we only process the first 0 ... i elements of this row,
   // as that's all that are kept in the symmetric packed format.
 
   Matrix<Real> temp_this(*this);
   Real *data = this->data_;
   const Real *Mdata = M.Data(), *MAdata = temp_MA.Data();
   MatrixIndexT temp_MA_stride = temp_MA.Stride(), Mstride = M.Stride();
 
   if (transM == kNoTrans) {
     // The column of M^T corresponds to the rows of the supplied matrix.
     for (MatrixIndexT i = 0; i < dim; i++, data += i) {
       MatrixIndexT num_rows = i + 1, num_cols = Adim;
       Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
                       temp_MA_stride, Mdata + (i * Mstride), 1, beta, data, 1);
     }
   } else {
     // The column of M^T corresponds to the columns of the supplied matrix.
     for (MatrixIndexT i = 0; i < dim; i++, data += i) {
       MatrixIndexT num_rows = i + 1, num_cols = Adim;
       Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
                       temp_MA_stride, Mdata + i, Mstride, beta, data, 1);
     }
   }
 }

◆ AddSp()

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< float >::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().

                                                          {
     this->AddPacked(alpha, Ma);
   }

◆ AddTp2()

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 SpMatrix< float >::AddSp(), CompressedAffineXformStats::ExtractOneG(), main(), and CuMatrixBase< float >::SymInvertPosDef().

                                                                          {
   Matrix<Real> Tmat(T);
   AddMat2(alpha, Tmat, transM, beta);
 }

◆ AddTp2Sp()

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 SpMatrix< float >::AddSp(), Sgmm2Project::ComputeLdaTransform(), kaldi::UnitTestPldaEstimation(), and PldaUnsupervisedAdaptor::UpdatePlda().

                                                {
   Matrix<Real> Tmat(T);
   AddMat2Sp(alpha, Tmat, transM, A, beta);
 }

◆ AddVec2() [1/4]

void AddVec2	(	const double	alpha,
		const VectorBase< double > &	v
	)

◆ AddVec2() [2/4]

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.

Referenced by CovarianceStats::AccStats(), IvectorExtractorUtteranceStats::AccStats(), LdaEstimate::Accumulate(), FmllrSgmm2Accs::AccumulateForFmllrSubspace(), RegtreeMllrDiagGmmAccs::AccumulateForGaussian(), RegtreeFmllrDiagGmmAccs::AccumulateForGaussian(), RegtreeMllrDiagGmmAccs::AccumulateForGmm(), RegtreeFmllrDiagGmmAccs::AccumulateForGmm(), AccumFullGmm::AccumulateFromPosteriors(), FmllrSgmm2Accs::AccumulateFromPosteriors(), SpMatrix< float >::AddDiagVec(), SpMatrix< float >::AddSp(), SpMatrix< float >::AddVec2(), FmllrRawAccs::CommitSingleFrameStats(), FmllrDiagGmmAccs::CommitSingleFrameStats(), IvectorExtractorStats::CommitStatsForM(), IvectorExtractorStats::CommitStatsForPrior(), IvectorExtractorStats::CommitStatsForWPoint(), Fmpe::ComputeC(), kaldi::ComputeFeatureNormalizingTransform(), AmSgmm2::ComputeFmllrPreXform(), Sgmm2Project::ComputeLdaStats(), EbwAmSgmm2Updater::ComputePhoneVecStats(), kaldi::EstPca(), 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().

                                                                              {
   KALDI_ASSERT(v.Dim() == this->NumRows());
   Real *data = this->data_;
   const OtherReal *v_data = v.Data();
   MatrixIndexT nr = this->num_rows_;
   for (MatrixIndexT i = 0; i < nr; i++)
     for (MatrixIndexT j = 0; j <= i; j++, data++)
       *data += alpha * v_data[i] * v_data[j];
 }

◆ AddVec2() [3/4]

void AddVec2	(	const float	alpha,
		const VectorBase< float > &	v
	)

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

                                                                            {
   KALDI_ASSERT(v.Dim() == this->NumRows());
   cblas_Xspr(v.Dim(), alpha, v.Data(), 1,
              this->data_);
 }

◆ AddVec2() [4/4]

void AddVec2	(	const double	alpha,
		const VectorBase< double > &	v
	)

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

                                                                               {
   KALDI_ASSERT(v.Dim() == num_rows_);
   cblas_Xspr(v.Dim(), alpha, v.Data(), 1, data_);
 }

◆ AddVec2Sp()

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.

Referenced by SpMatrix< float >::AddSp(), kaldi::SolveQuadraticMatrixProblem(), kaldi::SolveQuadraticProblem(), and kaldi::UnitTestAddVec2Sp().

                                                                          {
   KALDI_ASSERT(v.Dim() == this->NumRows() && S.NumRows() == this->NumRows());
   const Real *Sdata = S.Data();
   const Real *vdata = v.Data();
   Real *data = this->data_;
   MatrixIndexT dim = this->num_rows_;
   for (MatrixIndexT r = 0; r < dim; r++)
     for (MatrixIndexT c = 0; c <= r; c++, Sdata++, data++)
       *data = beta * *data + alpha * vdata[r] * vdata[c] * *Sdata;
 }

◆ AddVecVec()

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.

Referenced by SpMatrix< float >::AddSp(), and kaldi::UnitTestSpAddVecVec().

                                                           {
   int32 dim = this->NumRows();
   KALDI_ASSERT(dim == v.Dim() && dim == w.Dim() && dim > 0);
   cblas_Xspr2(dim, alpha, v.Data(), 1, w.Data(), 1, this->data_);
 }

◆ ApplyFloor() [1/2]

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.

Referenced by SpMatrix< float >::AddSp(), main(), kaldi::UnitTestFloorChol(), kaldi::UnitTestFloorUnit(), IvectorExtractorStats::UpdateVariances(), and EbwAmSgmm2Updater::UpdateVars().

                                              {
   MatrixIndexT dim = this->NumRows();
   int nfloored = 0;
   KALDI_ASSERT(C.NumRows() == dim);
   KALDI_ASSERT(alpha > 0);
   TpMatrix<Real> L(dim);
   L.Cholesky(C);
   L.Scale(std::sqrt(alpha));  // equivalent to scaling C by alpha.
   TpMatrix<Real> LInv(L);
   LInv.Invert();
 
   SpMatrix<Real> D(dim);
   {  // D = L^{-1} * (*this) * L^{-T}
     Matrix<Real> LInvFull(LInv);
     D.AddMat2Sp(1.0, LInvFull, kNoTrans, (*this), 0.0);
   }
 
   Vector<Real> l(dim);
   Matrix<Real> U(dim, dim);
 
   D.Eig(&l, &U);
 
   if (verbose) {
     KALDI_LOG << "ApplyFloor: flooring following diagonal to 1: " << l;
   }
   for (MatrixIndexT i = 0; i < l.Dim(); i++) {
     if (l(i) < 1.0) {
       nfloored++;
       l(i) = 1.0;
     }
   }
   l.ApplyPow(0.5);
   U.MulColsVec(l);
   D.AddMat2(1.0, U, kNoTrans, 0.0);
   {  // D' := U * diag(l') * U^T ... l'=floor(l, 1)
     Matrix<Real> LFull(L);
     (*this).AddMat2Sp(1.0, LFull, kNoTrans, D, 0.0);  // A := L * D' * L^T
   }
   return nfloored;
 }

◆ ApplyFloor() [2/2]

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.

                                          {
   MatrixIndexT Dim = this->NumRows();
   int nfloored = 0;
   Vector<Real> s(Dim);
   Matrix<Real> P(Dim, Dim);
   (*this).Eig(&s, &P);
   for (MatrixIndexT i = 0; i < Dim; i++) {
     if (s(i) < floor) {
       nfloored++;
       s(i) = floor;
     }
   }
   (*this).AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
   return nfloored;
 }

◆ ApplyPow()

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.

Referenced by SpMatrix< float >::Cond(), AmSgmm2::SplitSubstates(), and kaldi::UnitTestPower().

                                         {
   if (power == 1) return;  // can do nothing.
   MatrixIndexT D = this->NumRows();
   KALDI_ASSERT(D > 0);
   Matrix<Real> U(D, D);
   Vector<Real> l(D);
   (*this).SymPosSemiDefEig(&l, &U);
 
   Vector<Real> l_copy(l);
   try {
     l.ApplyPow(power * 0.5);
   }
   catch(...) {
     KALDI_ERR << "Error taking power " << (power * 0.5) << " of vector "
               << l_copy;
   }
   U.MulColsVec(l);
   (*this).AddMat2(1.0, U, kNoTrans, 0.0);
 }

◆ ApproxEqual()

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.

Referenced by SpMatrix< float >::AddSp(), kaldi::AssertEqual(), kaldi::UnitTestEigSp(), and kaldi::UnitTestNorm().

                                                                              {
   if (this->NumRows() != other.NumRows())
     KALDI_ERR << "SpMatrix::AproxEqual, size mismatch, "
               << this->NumRows() << " vs. " << other.NumRows();
   SpMatrix<Real> tmp(*this);
   tmp.AddSp(-1.0, other);
   return (tmp.FrobeniusNorm() <= tol * std::max(this->FrobeniusNorm(), other.FrobeniusNorm()));
 }

◆ Cond()

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

                            {
     Matrix<Real> tmp(*this);
     return tmp.Cond();
   }

◆ CopyFromMat()

void CopyFromMat	(	const MatrixBase< Real > &	orig,
		SpCopyType	copy_type = `kTakeMean`
	)

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

Referenced by BasisFmllrEstimate::ComputeAmDiagPrecond(), SpMatrix< float >::CopyFromSp(), SpMatrix< float >::SpMatrix(), kaldi::UnitTestCopySp(), kaldi::UnitTestDeterminant(), kaldi::UnitTestDeterminantSign(), and kaldi::UnitTestPower().

                                                        {
   KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
   MatrixIndexT D = this->NumRows();
 
   switch (copy_type) {
     case kTakeMeanAndCheck:
       {
         Real good_sum = 0.0, bad_sum = 0.0;
         for (MatrixIndexT i = 0; i < D; i++) {
           for (MatrixIndexT j = 0; j < i; j++) {
             Real a = M(i, j), b = M(j, i), avg = 0.5*(a+b), diff = 0.5*(a-b);
             (*this)(i, j) = avg;
             good_sum += std::abs(avg);
             bad_sum += std::abs(diff);
           }
           good_sum += std::abs(M(i, i));
           (*this)(i, i) = M(i, i);
         }
         if (bad_sum > 0.01 * good_sum) {
           KALDI_ERR << "SpMatrix::Copy(), source matrix is not symmetric: "
                     << bad_sum <<  ">" << good_sum;
         }
         break;
       }
     case kTakeMean:
       {
         for (MatrixIndexT i = 0; i < D; i++) {
           for (MatrixIndexT j = 0; j < i; j++) {
             (*this)(i, j) = 0.5*(M(i, j) + M(j, i));
           }
           (*this)(i, i) = M(i, i);
         }
         break;
       }
     case kTakeLower:
       { // making this one a bit more efficient.
         const Real *src = M.Data();
         Real *dest = this->data_;
         MatrixIndexT stride = M.Stride();
         for (MatrixIndexT i = 0; i < D; i++) {
           for (MatrixIndexT j = 0; j <= i; j++)
             dest[j] = src[j];
           dest += i + 1;
           src += stride;
         }
       }
       break;
     case kTakeUpper:
       for (MatrixIndexT i = 0; i < D; i++)
         for (MatrixIndexT j = 0; j <= i; j++)
           (*this)(i, j) = M(j, i);
       break;
     default:
       KALDI_ASSERT("Invalid argument to SpMatrix::CopyFromMat");
   }
 }

◆ CopyFromSp() [1/2]

void CopyFromSp ( const SpMatrix< Real > & other )

inline

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

Referenced by SpMatrix< float >::AddMat2Sp(), PldaEstimator::ComputeObjfPart2(), DiagGmm::CopyFromFullGmm(), PldaEstimator::EstimateFromStats(), FullGmm::GetComponentMean(), AmSgmm2::GetInvCovars(), IvectorExtractor::GetIvectorDistribution(), FullGmm::GetMeans(), LdaEstimate::GetStats(), PldaEstimator::GetStatsFromClassMeans(), main(), FullGmm::MergedComponentsLogdet(), Sgmm2Project::ProjectVariance(), AffineXformStats::Read(), FullGmm::SetInvCovars(), kaldi::UnitTestIoCross(), kaldi::UnitTestMat2Vec(), kaldi::UnitTestSpInvert(), and MleAmSgmm2Updater::UpdateVars().

                                                {
     PackedMatrix<Real>::CopyFromPacked(other);
   }

◆ CopyFromSp() [2/2]

void CopyFromSp ( const SpMatrix< OtherReal > & other )

inline

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

                                                     {
     PackedMatrix<Real>::CopyFromPacked(other);
   }

◆ Eig()

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.

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

                                                                        {
   MatrixIndexT dim = this->NumRows();
   KALDI_ASSERT(s->Dim() == dim);
   KALDI_ASSERT(P == NULL || (P->NumRows() == dim && P->NumCols() == dim));
 
   SpMatrix<Real> A(*this); // Copy *this, since the tridiagonalization
   // and QR decomposition are destructive.
   // Note: for efficiency of memory access, the tridiagonalization
   // algorithm makes the *rows* of P the eigenvectors, not the columns.
   // We'll transpose P before we exit.
   // Also note: P may be null if you don't want the eigenvectors.  This
   // will make this function more efficient.
 
   A.Tridiagonalize(P); // Tridiagonalizes.
   A.Qr(P); // Diagonalizes.
   if(P) P->Transpose();
   s->CopyDiagFromPacked(A);
 }

◆ EigInternal()

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

private

Referenced by SpMatrix< float >::LimitCondDouble().

◆ FrobeniusNorm()

Real FrobeniusNorm ( ) const

sqrt of sum of square elements.

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

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

                                          {
   Real sum = 0.0;
   MatrixIndexT R = this->NumRows();
   for (MatrixIndexT i = 0; i < R; i++) {
     for (MatrixIndexT j = 0; j < i; j++)
       sum += (*this)(i, j) * (*this)(i, j) * 2;
     sum += (*this)(i, i) * (*this)(i, i);
   }
   return std::sqrt(sum);
 }

◆ Invert()

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.

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

                                                                            {
   // these are CLAPACK types
   KaldiBlasInt   result;
   KaldiBlasInt   rows = static_cast<int>(this->num_rows_);
   KaldiBlasInt*  p_ipiv = new KaldiBlasInt[rows];
   Real *p_work;  // workspace for the lapack function
   void *temp;
   if ((p_work = static_cast<Real*>(
           KALDI_MEMALIGN(16, sizeof(Real) * rows, &temp))) == NULL) {
     delete[] p_ipiv;
     throw std::bad_alloc();
   }
 #ifdef HAVE_OPENBLAS
   memset(p_work, 0, sizeof(Real) * rows); // gets rid of a probably
   // spurious Valgrind warning about jumps depending upon uninitialized values.
 #endif
 
 
   // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
   // of the data. We have a row-wise storage, therefore, we need to "invert"
   clapack_Xsptrf(&rows, this->data_, p_ipiv, &result);
 
 
   KALDI_ASSERT(result >= 0 && "Call to CLAPACK ssptrf_ called with wrong arguments");
 
   if (result > 0) {  // Singular...
     if (det_sign) *det_sign = 0;
     if (logdet) *logdet = -std::numeric_limits<Real>::infinity();
     if (need_inverse) KALDI_ERR << "CLAPACK stptrf_ : factorization failed";
   } else {  // Not singular.. compute log-determinant if needed.
     if (logdet != NULL || det_sign != NULL) {
       Real prod = 1.0, log_prod = 0.0;
       int sign = 1;
       for (int i = 0; i < (int)this->num_rows_; i++) {
         if (p_ipiv[i] > 0) {  // not a 2x2 block...
           // if (p_ipiv[i] != i+1) sign *= -1;  // row swap.
           Real diag = (*this)(i, i);
           prod *= diag;
         } else {  // negative: 2x2 block. [we are in first of the two].
           i++;  // skip over the first of the pair.
           // each 2x2 block...
           Real diag1 = (*this)(i, i), diag2 = (*this)(i-1, i-1),
               offdiag = (*this)(i, i-1);
           Real thisdet = diag1*diag2 - offdiag*offdiag;
           // thisdet == determinant of 2x2 block.
           // The following line is more complex than it looks: there are 2 offsets of
           // 1 that cancel.
           prod *= thisdet;
         }
         if (i == (int)(this->num_rows_-1) || fabs(prod) < 1.0e-10 || fabs(prod) > 1.0e+10) {
           if (prod < 0) { prod = -prod; sign *= -1; }
           log_prod += kaldi::Log(std::abs(prod));
           prod = 1.0;
         }
       }
       if (logdet != NULL) *logdet = log_prod;
       if (det_sign != NULL) *det_sign = sign;
     }
   }
   if (!need_inverse) {
     delete [] p_ipiv;
     KALDI_MEMALIGN_FREE(p_work);
     return;  // Don't need what is computed next.
   }
   // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
   // of the data. We have a row-wise storage, therefore, we need to "invert"
   clapack_Xsptri(&rows, this->data_, p_ipiv, p_work, &result);
 
   KALDI_ASSERT(result >=0 &&
                "Call to CLAPACK ssptri_ called with wrong arguments");
 
   if (result != 0) {
     KALDI_ERR << "CLAPACK ssptrf_ : Matrix is singular";
   }
 
   delete [] p_ipiv;
   KALDI_MEMALIGN_FREE(p_work);
 }

◆ InvertDouble()

void InvertDouble	(	Real *	logdet = `NULL`,
		Real *	det_sign = `NULL`,
		bool	inverse_needed = `true`
	)

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

Referenced by AmSgmm2::ComputeFmllrPreXform(), FullGmm::ComputeGconsts(), FullGmm::GetComponentMean(), FullGmm::GetMeans(), main(), SpMatrix< float >::operator=(), FullGmm::SetInvCovars(), and kaldi::UnitTestSpInvert().

                                                        {
   SpMatrix<double> dmat(*this);
   double logdet_tmp, det_sign_tmp;
   dmat.Invert(logdet ? &logdet_tmp : NULL,
               det_sign ? &det_sign_tmp : NULL,
               inverse_needed);
   if (logdet) *logdet = logdet_tmp;
   if (det_sign) *det_sign = det_sign_tmp;
   (*this).CopyFromSp(dmat);
 }

◆ IsDiagonal()

bool IsDiagonal ( Real cutoff = 1.0e-05 ) const

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

Referenced by SpMatrix< float >::AddSp(), AmSgmm2::ComputeFmllrPreXform(), MleAmSgmm2Updater::RenormalizeV(), kaldi::SolveDoubleQuadraticMatrixProblem(), kaldi::UnitTestPca2(), and kaldi::UnitTestTridiagonalizeAndQr().

                                                  {
   MatrixIndexT R = this->NumRows();
   Real bad_sum = 0.0, good_sum = 0.0;
   for (MatrixIndexT i = 0; i < R; i++) {
     for (MatrixIndexT j = 0; j <= i; j++) {
       if (i == j)
         good_sum += std::abs((*this)(i, j));
       else
         bad_sum += std::abs((*this)(i, j));
     }
   }
   return (!(bad_sum > good_sum * cutoff));
 }

◆ IsPosDef()

bool IsPosDef ( ) const

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

Referenced by SpMatrix< float >::PrintEigs(), and MleAmSgmm2Updater::RenormalizeV().

                                     {
   MatrixIndexT D = (*this).NumRows();
   KALDI_ASSERT(D > 0);
   try {
     TpMatrix<Real> C(D);
     C.Cholesky(*this);
     for (MatrixIndexT r = 0; r < D; r++)
       if (C(r, r) == 0.0) return false;
     return true;
   }
   catch(...) {  // not positive semidefinite.
     return false;
   }
 }

◆ IsTridiagonal()

bool IsTridiagonal ( Real cutoff = 1.0e-05 ) const

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

Referenced by SpMatrix< float >::AddSp(), and kaldi::UnitTestTridiag().

                                                     {
   MatrixIndexT R = this->NumRows();
   Real max_abs_2diag = 0.0, max_abs_offdiag = 0.0;
   for (MatrixIndexT i = 0; i < R; i++)
     for (MatrixIndexT j = 0; j <= i; j++) {
       if (j+1 < i)
         max_abs_offdiag = std::max(max_abs_offdiag,
                                    std::abs((*this)(i, j)));
       else
         max_abs_2diag = std::max(max_abs_2diag,
                                  std::abs((*this)(i, j)));
     }
   return (max_abs_offdiag <= cutoff * max_abs_2diag);
 }

◆ IsUnit()

bool IsUnit ( Real cutoff = 1.0e-05 ) const

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

Referenced by SpMatrix< float >::AddSp(), AmSgmm2::ComputeFmllrPreXform(), Sgmm2Project::ComputeLdaTransform(), PldaEstimator::GetOutput(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), MleAmSgmm2Updater::RenormalizeV(), OnlineNaturalGradient::ReorthogonalizeRt1(), OnlinePreconditioner::ReorthogonalizeXt1(), OnlinePreconditioner::SelfTest(), OnlineNaturalGradient::SelfTest(), kaldi::SolveDoubleQuadraticMatrixProblem(), kaldi::UnitTestCuSpMatrixConstructor(), kaldi::UnitTestPca(), kaldi::UnitTestTopEigs(), and IvectorExtractorStats::UpdatePrior().

                                              {
   MatrixIndexT R = this->NumRows();
   Real max = 0.0;  // max error
   for (MatrixIndexT i = 0; i < R; i++)
     for (MatrixIndexT j = 0; j <= i; j++)
       max = std::max(max, static_cast<Real>(std::abs((*this)(i, j) -
                                                      (i == j ? 1.0 : 0.0))));
   return (max <= cutoff);
 }

◆ IsZero()

bool IsZero ( Real cutoff = 1.0e-05 ) const

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

Referenced by SpMatrix< float >::AddSp(), kaldi::SolveQuadraticMatrixProblem(), kaldi::SolveQuadraticProblem(), and kaldi::UnitTestResize().

                                              {
   if (this->num_rows_ == 0) return true;
   return (this->Max() <= cutoff && this->Min() >= -cutoff);
 }

◆ LimitCond()

MatrixIndexT LimitCond	(	Real	maxCond = `1.0e+5`,
		bool	invert = `false`
	)

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

Referenced by SpMatrix< float >::AddSp(), SpMatrix< float >::LimitCondDouble(), kaldi::UnitTestLimitCond(), and kaldi::UnitTestLimitCondInvert().

                                                                 {  // e.g. maxCond = 1.0e+05.
   MatrixIndexT Dim = this->NumRows();
   Vector<Real> s(Dim);
   Matrix<Real> P(Dim, Dim);
   (*this).SymPosSemiDefEig(&s, &P);
   KALDI_ASSERT(maxCond > 1);
   Real floor = s.Max() / maxCond;
   if (floor < 0) floor = 0;
   if (floor < 1.0e-40) {
     KALDI_WARN << "LimitCond: limiting " << floor << " to 1.0e-40";
     floor = 1.0e-40;
   }
   MatrixIndexT nfloored = 0;
   for (MatrixIndexT i = 0; i < Dim; i++) {
     if (s(i) <= floor) nfloored++;
     if (invert)
       s(i) = 1.0 / std::sqrt(std::max(s(i), floor));
     else
       s(i) = std::sqrt(std::max(s(i), floor));
   }
   P.MulColsVec(s);
   (*this).AddMat2(1.0, P, kNoTrans, 0.0);  // (*this) = P*P^T.  ... (*this) = P * floor(s) * P^T  ... if P was original P.
   return nfloored;
 }

◆ LimitCondDouble()

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

                                                                            {
     SpMatrix<double> dmat(*this);
     MatrixIndexT ans = dmat.LimitCond(maxCond, invert);
     (*this).CopyFromSp(dmat);
     return ans;
   }

◆ LogDet()

Real LogDet ( Real * det_sign = NULL ) const

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

Referenced by SpMatrix< float >::AddSp(), kaldi::UnitTestDeterminant(), kaldi::UnitTestDeterminantSign(), kaldi::UnitTestTridiagonalize(), kaldi::UnitTestTridiagonalizeAndQr(), and EbwAmSgmm2Updater::UpdateVars().

                                                 {
   Real log_det;
   SpMatrix<Real> tmp(*this);
   // false== output not needed (saves some computation).
   tmp.Invert(&log_det, det_sign, false);
   return log_det;
 }

◆ LogPosDefDet()

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.

Referenced by SpMatrix< float >::AddSp(), FullGmm::ComputeGconsts(), kaldi::GetLogDetNoFailure(), FullGmm::MergedComponentsLogdet(), IvectorExtractorStats::PriorDiagnostics(), rand_posdef_spmatrix(), kaldi::unittest::RandPosdefSpMatrix(), RandPosdefSpMatrix(), kaldi::UnitTestDeterminant(), kaldi::UnitTestDeterminantSign(), and IvectorExtractorStats::UpdateVariances().

                                         {
   TpMatrix<Real> chol(this->NumRows());
   double det = 0.0;
   double diag;
   chol.Cholesky(*this);  // Will throw exception if not +ve definite!
 
   for (MatrixIndexT i = 0; i < this->NumRows(); i++) {
     diag = static_cast<double>(chol(i, i));
     det += kaldi::Log(diag);
   }
   return static_cast<Real>(2*det);
 }

◆ MaxAbsEig()

Real MaxAbsEig ( ) const

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

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

Referenced by SpMatrix< float >::Cond(), main(), kaldi::UnitTestLinearCgd(), kaldi::UnitTestMaxAbsEig(), and IvectorExtractorStats::UpdateVariances().

                                      {
   Vector<Real> s(this->NumRows());
   this->Eig(&s, static_cast<MatrixBase<Real>*>(NULL));
   return std::max(s.Max(), -s.Min());
 }

◆ operator()() [1/2]

Real operator()	(	MatrixIndexT	r,
		MatrixIndexT	c
	)		const

inline

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

                                                                 {
     // if column is less than row, then swap these as matrix is stored
     // as upper-triangular...  only allowed for const matrix object.
     if (static_cast<UnsignedMatrixIndexT>(c) >
         static_cast<UnsignedMatrixIndexT>(r))
       std::swap(c, r);
     // c<=r now so don't have to check c.
     KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
                  static_cast<UnsignedMatrixIndexT>(this->num_rows_));
     return *(this->data_ + (r*(r+1)) / 2 + c);
     // Duplicating code from PackedMatrix.h
   }

◆ operator()() [2/2]

Real& operator()	(	MatrixIndexT	r,
		MatrixIndexT	c
	)

inline

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

                                                            {
     if (static_cast<UnsignedMatrixIndexT>(c) >
         static_cast<UnsignedMatrixIndexT>(r))
       std::swap(c, r);
     // c<=r now so don't have to check c.
     KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
                  static_cast<UnsignedMatrixIndexT>(this->num_rows_));
     return *(this->data_ + (r * (r + 1)) / 2 + c);
     // Duplicating code from PackedMatrix.h
   }

◆ operator=()

SpMatrix<Real>& operator= ( const SpMatrix< Real > & other )

inline

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

                                                          {
     PackedMatrix<Real>::operator=(other);
     return *this;
   }

◆ PrintEigs()

void PrintEigs ( const char * name )

inline

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

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

                                    {
     Vector<Real> s((*this).NumRows());
     Matrix<Real> P((*this).NumRows(), (*this).NumCols());
     SymPosSemiDefEig(&s, &P);
     KALDI_LOG << "PrintEigs: " << name << ": " << s;
   }

◆ Qr()

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.

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

                                            {
   KALDI_ASSERT(this->IsTridiagonal());
   // We envisage that Q would be square but we don't check for this,
   // as there are situations where you might not want this.
   KALDI_ASSERT(Q == NULL || Q->NumRows() == this->NumRows());
   // Note: the first couple of lines of the algorithm they give would be done
   // outside of this function, by calling Tridiagonalize().
 
   MatrixIndexT n = this->NumRows();
   Vector<Real> diag(n), off_diag(n-1);
   for (MatrixIndexT i = 0; i < n; i++) {
     diag(i) = (*this)(i, i);
     if (i > 0) off_diag(i-1) = (*this)(i, i-1);
   }
   QrInternal(n, diag.Data(), off_diag.Data(), Q);
   // Now set *this to the value represented by diag and off_diag.
   this->SetZero();
   for (MatrixIndexT i = 0; i < n; i++) {
     (*this)(i, i) = diag(i);
     if (i > 0) (*this)(i, i-1) = off_diag(i-1);
   }
 }

◆ Resize()

void Resize	(	MatrixIndexT	nRows,
		MatrixResizeType	resize_type = `kSetZero`
	)

inline

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

Referenced by SpMatrix< float >::AddMat2Sp(), BasisFmllrEstimate::ComputeAmDiagPrecond(), AmSgmm2::ComputeHsmFromModel(), Sgmm2Project::ComputeLdaStats(), FastNnetCombiner::ComputePreconditioner(), Sgmm2Project::ComputeProjection(), FmllrRawAccs::ConvertToSimpleStats(), kaldi::EstPca(), FisherComputationClass::FisherComputationClass(), FmllrRawAccs::FmllrRawAccs(), AmSgmm2::GetInvCovars(), LdaEstimate::GetStats(), LdaEstimate::Init(), PldaEstimator::InitParameters(), IvectorExtractorStats::IvectorExtractorStats(), main(), Sgmm2Project::ProjectVariance(), AffineXformStats::Read(), PldaEstimator::ResetPerIterStats(), kaldi::UnitTestIo(), kaldi::UnitTestResize(), and FisherComputationClass::~FisherComputationClass().

                                                                                   {
     PackedMatrix<Real>::Resize(nRows, resize_type);
   }

◆ Swap()

void Swap ( SpMatrix< Real > * other )

Shallow swap.

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

Referenced by SpMatrix< float >::SpMatrix().

                                                {
   std::swap(this->data_, other->data_);
   std::swap(this->num_rows_, other->num_rows_);
 }

◆ SymPosSemiDefEig()

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.

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

                                                             {
   Eig(s, P);
   Real max = s->Max(), min = s->Min();
   KALDI_ASSERT(-min <= tolerance * max);
   s->ApplyFloor(0.0);
 }

◆ TopEigs()

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.

Referenced by kaldi::ComputePca(), SpMatrix< float >::Cond(), SpMatrix< float >::TopEigs(), and kaldi::UnitTestTopEigs().

                                                              {
   const SpMatrix<Real> &S(*this); // call this "S" for easy notation.
   MatrixIndexT eig_dim = s->Dim(); // Space of dim we want to retain.
   if (lanczos_dim <= 0)
     lanczos_dim = std::max(eig_dim + 50, eig_dim + eig_dim/2);
   MatrixIndexT dim = this->NumRows();
   if (lanczos_dim >= dim) {
     // There would be no speed advantage in using this method, so just
     // use the regular approach.
     Vector<Real> s_tmp(dim);
     Matrix<Real> P_tmp(dim, dim);
     this->Eig(&s_tmp, &P_tmp);
     SortSvd(&s_tmp, &P_tmp);
     s->CopyFromVec(s_tmp.Range(0, eig_dim));
     P->CopyFromMat(P_tmp.Range(0, dim, 0, eig_dim));
     return;
   }
   KALDI_ASSERT(eig_dim <= dim && eig_dim > 0);
   KALDI_ASSERT(P->NumRows() == dim && P->NumCols() == eig_dim); // each column
   // is one eigenvector.
 
   Matrix<Real> Q(lanczos_dim, dim); // The rows of Q will be the
   // orthogonal vectors of the Krylov subspace.
 
   SpMatrix<Real> T(lanczos_dim); // This will be equal to Q S Q^T,
   // i.e. *this projected into the Krylov subspace.  Note: only the
   // diagonal and off-diagonal fo T are nonzero, i.e. it's tridiagonal,
   // but we don't have access to the low-level algorithms that work
   // on that type of matrix (since we want to use ATLAS).  So we just
   // do normal SVD, on a full matrix; it won't typically dominate.
 
   Q.Row(0).SetRandn();
   Q.Row(0).Scale(1.0 / Q.Row(0).Norm(2));
   for (MatrixIndexT d = 0; d < lanczos_dim; d++) {
     Vector<Real> r(dim);
     r.AddSpVec(1.0, S, Q.Row(d), 0.0);
     // r = S * q_d
     MatrixIndexT counter = 0;
     Real end_prod;
     while (1) { // Normally we'll do this loop only once:
       // we repeat to handle cases where r gets very much smaller
       // and we want to orthogonalize again.
       // We do "full orthogonalization" to preserve stability,
       // even though this is usually a waste of time.
       Real start_prod = VecVec(r, r);
       for (SignedMatrixIndexT e = d; e >= 0; e--) { // e must be signed!
         SubVector<Real> q_e(Q, e);
         Real prod = VecVec(r, q_e);
         if (counter == 0 && static_cast<MatrixIndexT>(e) + 1 >= d) // Keep T tridiagonal, which
           T(d, e) = prod; // mathematically speaking, it is.
         r.AddVec(-prod, q_e); // Subtract component in q_e.
       }
       if (d+1 == lanczos_dim) break;
       end_prod = VecVec(r, r);
       if (end_prod <= 0.1 * start_prod) {
         // also handles case where both are 0.
         // We're not confident any more that it's completely
         // orthogonal to the rest so we want to re-do.
         if (end_prod == 0.0)
           r.SetRandn(); // "Restarting".
         counter++;
         if (counter > 100)
           KALDI_ERR << "Loop detected in Lanczos iteration.";
       } else {
         break;
       }
     }
     if (d+1 != lanczos_dim) {
       // OK, at this point we're satisfied that r is orthogonal
       // to all previous rows.
       KALDI_ASSERT(end_prod != 0.0); // should have looped.
       r.Scale(1.0 / std::sqrt(end_prod)); // make it unit.
       Q.Row(d+1).CopyFromVec(r);
     }
   }
 
   Matrix<Real> R(lanczos_dim, lanczos_dim);  
   R.SetUnit();
   T.Qr(&R); // Diagonalizes T.
   Vector<Real> s_tmp(lanczos_dim);
   s_tmp.CopyDiagFromSp(T);  
 
   // Now T = R * diag(s_tmp) * R^T.
   // The next call sorts the elements of s from greatest to least absolute value,
   // and moves around the rows of R in the corresponding way.  This picks out
   // the largest (absolute) eigenvalues.
   SortSvd(&s_tmp, static_cast<Matrix<Real>*>(NULL), &R);
   // Keep only the initial rows of R, those corresponding to greatest (absolute)
   // eigenvalues.
   SubMatrix<Real> Rsub(R, 0, eig_dim, 0, lanczos_dim);
   SubVector<Real> s_sub(s_tmp, 0, eig_dim);
   s->CopyFromVec(s_sub);
       
   // For working out what to do now, just assume the other eigenvalues were
   // zero.  This is just for purposes of knowing how to get the result, and
   // not getting things wrongly transposed.
   // We have T = Rsub^T * diag(s_sub) * Rsub.
   // Now, T = Q S Q^T, with Q orthogonal,  so S = Q^T T Q = Q^T Rsub^T * diag(s) * Rsub * Q.
   // The output is P and we want S = P * diag(s) * P^T, so we need P = Q^T Rsub^T.
   P->AddMatMat(1.0, Q, kTrans, Rsub, kTrans, 0.0);
 }

◆ Trace()

Real Trace ( ) const

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

Referenced by FastNnetCombiner::ComputePreconditioner(), LdaEstimate::Estimate(), PldaEstimator::EstimateFromStats(), FeatureTransformEstimate::EstimateInternal(), IvectorExtractor::GetIvectorDistribution(), IvectorExtractor::GetPriorAuxf(), SpMatrix< float >::LimitCondDouble(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), OnlinePreconditioner::PreconditionDirectionsInternal(), OnlineNaturalGradient::PreconditionDirectionsInternal(), CompressedAffineXformStats::PrepareOneG(), IvectorExtractorStats::PriorDiagnostics(), kaldi::UnitTestCuMatrixSymInvertPosDef(), kaldi::UnitTestCuSpMatrixCopyFromMat(), kaldi::UnitTestPca2(), kaldi::UnitTestPldaEstimation(), kaldi::UnitTestTrace(), kaldi::UnitTestTridiagonalize(), and kaldi::UnitTestTridiagonalizeAndQr().

                                  {
   const Real *data = this->data_;
   MatrixIndexT num_rows = this->num_rows_;
   Real ans = 0.0;
   for (int32 i = 1; i <= num_rows; i++, data += i)
     ans += *data;
   return ans;
 }

◆ Tridiagonalize()

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.

Referenced by SpMatrix< float >::Eig(), SpMatrix< float >::LimitCondDouble(), SpMatrix< float >::Tridiagonalize(), kaldi::UnitTestTridiagonalize(), and kaldi::UnitTestTridiagonalizeAndQr().

                                                        {
   MatrixIndexT n = this->NumRows();
   KALDI_ASSERT(Q == NULL || (Q->NumRows() == n &&
                              Q->NumCols() == n));
   if (Q != NULL) Q->SetUnit();
   Real *data = this->Data();
   Real *qdata = (Q == NULL ? NULL : Q->Data());
   MatrixIndexT qstride = (Q == NULL ? 0 : Q->Stride());
   Vector<Real> tmp_v(n-1), tmp_p(n);
   Real beta, *v = tmp_v.Data(), *p = tmp_p.Data(), *w = p, *x = p;
   for (MatrixIndexT k = n-1; k >= 2; k--) {
     MatrixIndexT ksize = ((k+1)*k)/2;
     // ksize is the packed size of the lower-triangular matrix of size k,
     // which is the size of "all rows previous to this one."
     Real *Arow = data + ksize; // In Golub+Van Loan it was A(k+1:n, k), we
     // have Arow = A(k, 0:k-1).
     HouseBackward(k, Arow, v, &beta); // sets v and beta.
     cblas_Xspmv(k, beta, data, v, 1, 0.0, p, 1); // p = beta * A(0:k-1,0:k-1) v
     Real minus_half_beta_pv = -0.5 * beta * cblas_Xdot(k, p, 1, v, 1);
     cblas_Xaxpy(k, minus_half_beta_pv, v, 1, w, 1); // w = p - (beta p^T v/2) v;
     // this relies on the fact that w and p are the same pointer.
     // We're doing A(k, k-1) = ||Arow||.  It happens that this element
     // is indexed at ksize + k - 1 in the packed lower-triangular format.
     data[ksize + k - 1] = std::sqrt(cblas_Xdot(k, Arow, 1, Arow, 1));
     for (MatrixIndexT i = 0; i + 1 < k; i++)
       data[ksize + i] = 0; // This is not in Golub and Van Loan but is
     // necessary if we're not using parts of A to store the Householder
     // vectors.
     // We're doing A(0:k-1,0:k-1) -= (v w' + w v')
     cblas_Xspr2(k, -1.0, v, 1, w, 1, data);
     if (Q != NULL) { // C.f. Golub, Q is H_1 .. H_n-2... in this
       // case we apply them in the opposite order so it's H_n-1 .. H_1,
       // but also Q is transposed so we really have Q = H_1 .. H_n-1.
       // It's a double negative.    
       // Anyway, we left-multiply Q by each one.  The H_n would each be
       // diag(I + beta v v', I) but we don't ever touch the last dims.
       // We do (in Matlab notation):
       // Q(0:k-1,:) = (I - beta v v') * Q, i.e.:
       // Q(:,0:i-1) += -beta v (v' Q(:,0:k-1)v .. let x = -beta Q(0:k-1,:)^T v.
       cblas_Xgemv(kTrans, k, n, -beta, qdata, qstride, v, 1, 0.0, x, 1);
       // now x = -beta Q(:,0:k-1) v.
       // The next line does: Q(:,0:k-1) += v x'.
       cblas_Xger(k, n, 1.0, v, 1, x, 1, qdata, qstride);
     }
   }
 }

Friends And Related Function Documentation

◆ CuSpMatrix< Real >

friend class CuSpMatrix< Real >

friend

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

◆ std::vector< Matrix< Real > >

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:

matrix/matrix-common.h
matrix/sp-matrix.h
matrix/qr.cc
matrix/sp-matrix.cc
cudamatrix/cu-sp-matrix.h

Public Member Functions

Private Member Functions

Friends

Additional Inherited Members

Detailed Description

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

Constructor & Destructor Documentation

◆ SpMatrix() [1/6]

◆ SpMatrix() [2/6]

◆ SpMatrix() [3/6]

◆ SpMatrix() [4/6]

◆ SpMatrix() [5/6]

◆ SpMatrix() [6/6]

Member Function Documentation

◆ AddDiagVec()

◆ AddMat2()

◆ AddMat2Sp()

◆ AddMat2Vec()

◆ AddSmat2Sp()

◆ AddSp()

◆ AddTp2()

◆ AddTp2Sp()

◆ AddVec2() [1/4]

◆ AddVec2() [2/4]

◆ AddVec2() [3/4]

◆ AddVec2() [4/4]

◆ AddVec2Sp()

◆ AddVecVec()

◆ ApplyFloor() [1/2]

◆ ApplyFloor() [2/2]

◆ ApplyPow()

◆ ApproxEqual()

◆ Cond()

◆ CopyFromMat()

◆ CopyFromSp() [1/2]

◆ CopyFromSp() [2/2]

◆ Eig()

◆ EigInternal()

◆ FrobeniusNorm()

◆ Invert()

◆ InvertDouble()

◆ IsDiagonal()

◆ IsPosDef()

◆ IsTridiagonal()

◆ IsUnit()

◆ IsZero()

◆ LimitCond()

◆ LimitCondDouble()

◆ LogDet()

◆ LogPosDefDet()

◆ MaxAbsEig()

◆ operator()() [1/2]

◆ operator()() [2/2]

◆ operator=()

◆ PrintEigs()

◆ Qr()

◆ Resize()

◆ Swap()

◆ SymPosSemiDefEig()

◆ TopEigs()

◆ Trace()

◆ Tridiagonalize()

Friends And Related Function Documentation

◆ CuSpMatrix< Real >

◆ std::vector< Matrix< Real > >

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