All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Miscellaneous matrix/vector library functions and classes
Collaboration diagram for Miscellaneous matrix/vector library functions and classes:

Classes

class  MatrixExponential< Real >
 
struct  SolverOptions
 This class describes the options for maximizing various quadratic objective functions. More...
 
class  SplitRadixComplexFft< Real >
 
class  SplitRadixRealFft< Real >
 

Functions

template<typename Real >
void SortSvd (VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt=NULL, bool sort_on_absolute_value=true)
 Function to ensure that SVD is sorted. More...
 
template<typename Real >
void CreateEigenvalueMatrix (const VectorBase< Real > &real, const VectorBase< Real > &imag, MatrixBase< Real > *D)
 Creates the eigenvalue matrix D that is part of the decomposition used Matrix::Eig. More...
 
template<typename Real >
bool AttemptComplexPower (Real *x_re, Real *x_im, Real power)
 The following function is used in Matrix::Power, and separately tested, so we declare it here mainly for the testing code to see. More...
 
template<typename Real >
void ComplexFft (VectorBase< Real > *v, bool forward, Vector< Real > *tmp_work=NULL)
 The function ComplexFft does an Fft on the vector argument v. More...
 
template<typename Real >
void ComplexFt (const VectorBase< Real > &in, VectorBase< Real > *out, bool forward)
 ComplexFt is the same as ComplexFft but it implements the Fourier transform in an inefficient way. More...
 
template<typename Real >
void RealFft (VectorBase< Real > *v, bool forward)
 RealFft is a fourier transform of real inputs. More...
 
template<typename Real >
void RealFftInefficient (VectorBase< Real > *v, bool forward)
 Inefficient version of Fourier transform, for testing purposes. More...
 
template<typename Real >
void ComputeDctMatrix (Matrix< Real > *M)
 ComputeDctMatrix computes a matrix corresponding to the DCT, such that M * v equals the DCT of vector v. More...
 
template<typename Real >
void ComplexMul (const Real &a_re, const Real &a_im, Real *b_re, Real *b_im)
 ComplexMul implements, inline, the complex multiplication b *= a. More...
 
template<typename Real >
void ComplexAddProduct (const Real &a_re, const Real &a_im, const Real &b_re, const Real &b_im, Real *c_re, Real *c_im)
 ComplexMul implements, inline, the complex operation c += (a * b). More...
 
template<typename Real >
void ComplexImExp (Real x, Real *a_re, Real *a_im)
 ComplexImExp implements a <– exp(i x), inline. More...
 
template<typename Real >
void ComputePca (const MatrixBase< Real > &X, MatrixBase< Real > *U, MatrixBase< Real > *A, bool print_eigs=false, bool exact=true)
 ComputePCA does a PCA computation, using either outer products or inner products, whichever is more efficient. More...
 
template<typename Real >
void AddOuterProductPlusMinus (Real alpha, const VectorBase< Real > &a, const VectorBase< Real > &b, MatrixBase< Real > *plus, MatrixBase< Real > *minus)
 
template<typename Real1 , typename Real2 >
void AssertSameDim (const MatrixBase< Real1 > &mat1, const MatrixBase< Real2 > &mat2)
 
template<typename Real >
Real SolveQuadraticProblem (const SpMatrix< Real > &H, const VectorBase< Real > &g, const SolverOptions &opts, VectorBase< Real > *x)
 Maximizes the auxiliary function

\[ Q(x) = x.g - 0.5 x^T H x \]

using a numerically stable method. More...

 
template<typename Real >
Real SolveQuadraticMatrixProblem (const SpMatrix< Real > &Q, const MatrixBase< Real > &Y, const SpMatrix< Real > &P, const SolverOptions &opts, MatrixBase< Real > *M)
 Maximizes the auxiliary function :

\[ Q(x) = tr(M^T P Y) - 0.5 tr(P M Q M^T) \]

Like a numerically stable version of $ M := Y Q^{-1} $. More...

 
template<typename Real >
Real SolveDoubleQuadraticMatrixProblem (const MatrixBase< Real > &G, const SpMatrix< Real > &P1, const SpMatrix< Real > &P2, const SpMatrix< Real > &Q1, const SpMatrix< Real > &Q2, const SolverOptions &opts, MatrixBase< Real > *M)
 Maximizes the auxiliary function :

\[ Q(M) = tr(M^T G) -0.5 tr(P_1 M Q_1 M^T) -0.5 tr(P_2 M Q_2 M^T). \]

Encountered in matrix update with a prior. More...

 

Detailed Description

Function Documentation

void AddOuterProductPlusMinus ( Real  alpha,
const VectorBase< Real > &  a,
const VectorBase< Real > &  b,
MatrixBase< Real > *  plus,
MatrixBase< Real > *  minus 
)

Definition at line 933 of file matrix-functions.cc.

References VectorBase< Real >::Data(), MatrixBase< Real >::Data(), VectorBase< Real >::Dim(), rnnlm::i, rnnlm::j, KALDI_ASSERT, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), and MatrixBase< Real >::Stride().

Referenced by Fmpe::ApplyProjectionReverse(), and kaldi::UnitTestAddOuterProductPlusMinus().

937  {
938  KALDI_ASSERT(a.Dim() == plus->NumRows() && b.Dim() == plus->NumCols()
939  && a.Dim() == minus->NumRows() && b.Dim() == minus->NumCols());
940  int32 nrows = a.Dim(), ncols = b.Dim(), pskip = plus->Stride() - ncols,
941  mskip = minus->Stride() - ncols;
942  const Real *adata = a.Data(), *bdata = b.Data();
943  Real *plusdata = plus->Data(), *minusdata = minus->Data();
944 
945  for (int32 i = 0; i < nrows; i++) {
946  const Real *btmp = bdata;
947  Real multiple = alpha * *adata;
948  if (multiple > 0.0) {
949  for (int32 j = 0; j < ncols; j++, plusdata++, minusdata++, btmp++) {
950  if (*btmp > 0.0) *plusdata += multiple * *btmp;
951  else *minusdata -= multiple * *btmp;
952  }
953  } else {
954  for (int32 j = 0; j < ncols; j++, plusdata++, minusdata++, btmp++) {
955  if (*btmp < 0.0) *plusdata += multiple * *btmp;
956  else *minusdata -= multiple * *btmp;
957  }
958  }
959  plusdata += pskip;
960  minusdata += mskip;
961  adata++;
962  }
963 }
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void kaldi::AssertSameDim ( const MatrixBase< Real1 > &  mat1,
const MatrixBase< Real2 > &  mat2 
)
inline

Definition at line 222 of file matrix-functions.h.

References KALDI_ASSERT, MatrixBase< Real >::NumCols(), and MatrixBase< Real >::NumRows().

222  {
223  KALDI_ASSERT(mat1.NumRows() == mat2.NumRows()
224  && mat1.NumCols() == mat2.NumCols());
225 }
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
bool AttemptComplexPower ( Real *  x_re,
Real *  x_im,
Real  power 
)

The following function is used in Matrix::Power, and separately tested, so we declare it here mainly for the testing code to see.

It takes a complex value to a power using a method that will work for noninteger powers (but will fail if the complex value is real and negative).

Definition at line 2447 of file kaldi-matrix.cc.

Referenced by MatrixBase< Real >::Power(), and kaldi::UnitTestComplexPower().

2447  {
2448  // Used in Matrix<Real>::Power().
2449  // Attempts to take the complex value x to the power "power",
2450  // assuming that power is fractional (i.e. we don't treat integers as a
2451  // special case). Returns false if this is not possible, either
2452  // because x is negative and real (hence there is no obvious answer
2453  // that is "closest to 1", and anyway this case does not make sense
2454  // in the Matrix<Real>::Power() routine);
2455  // or because power is negative, and x is zero.
2456 
2457  // First solve for r and theta in
2458  // x_re = r*cos(theta), x_im = r*sin(theta)
2459  if (*x_re < 0.0 && *x_im == 0.0) return false; // can't do
2460  // it for negative real values.
2461  Real r = std::sqrt((*x_re * *x_re) + (*x_im * *x_im)); // r == radius.
2462  if (power < 0.0 && r == 0.0) return false;
2463  Real theta = std::atan2(*x_im, *x_re);
2464  // Take the power.
2465  r = std::pow(r, power);
2466  theta *= power;
2467  *x_re = r * std::cos(theta);
2468  *x_im = r * std::sin(theta);
2469  return true;
2470 }
void ComplexAddProduct ( const Real &  a_re,
const Real &  a_im,
const Real &  b_re,
const Real &  b_im,
Real *  c_re,
Real *  c_im 
)
inline

ComplexMul implements, inline, the complex operation c += (a * b).

Definition at line 38 of file matrix-functions-inl.h.

Referenced by kaldi::ComplexFftRecursive(), kaldi::ComplexFt(), SplitRadixRealFft< Real >::Compute(), and kaldi::RealFft().

40  {
41  *c_re += b_re*a_re - b_im*a_im;
42  *c_im += b_re*a_im + b_im*a_re;
43 }
void ComplexFft ( VectorBase< Real > *  v,
bool  forward,
Vector< Real > *  tmp_work = NULL 
)

The function ComplexFft does an Fft on the vector argument v.

v is a vector of even dimension, interpreted for both input and output as a vector of complex numbers i.e.

\[ v = ( re_0, im_0, re_1, im_1, ... ) \]

The dimension of v must be a power of 2.

If "forward == true" this routine does the Discrete Fourier Transform (DFT), i.e.:

\[ vout[m] \leftarrow \sum_{n = 0}^{N-1} vin[i] exp( -2pi m n / N ) \]

If "backward" it does the Inverse Discrete Fourier Transform (IDFT) WITHOUT THE FACTOR 1/N*, i.e.:

\[ vout[m] <-- \sum_{n = 0}^{N-1} vin[i] exp( 2pi m n / N ) \]

[note the sign difference on the 2 pi for the backward one.]

Note that this is the definition of the FT given in most texts, but it differs from the Numerical Recipes version in which the forward and backward algorithms are flipped.

Note that you would have to multiply by 1/N after the IDFT to get back to where you started from. We don't do this because in some contexts, the transform is made symmetric by multiplying by sqrt(N) in both passes. The user can do this by themselves.

See also SplitRadixComplexFft, declared in srfft.h, which is more efficient but only works if the length of the input is a power of 2.

Definition at line 334 of file matrix-functions.cc.

References kaldi::ComplexFftRecursive(), VectorBase< Real >::Data(), VectorBase< Real >::Dim(), kaldi::Factorize(), and KALDI_ASSERT.

Referenced by kaldi::RealFft(), kaldi::RealFftInefficient(), kaldi::UnitTestComplexFft(), and kaldi::UnitTestComplexFft2().

334  {
335  KALDI_ASSERT(v != NULL);
336 
337  if (v->Dim()<=1) return;
338  KALDI_ASSERT(v->Dim() % 2 == 0); // complex input.
339  int N = v->Dim() / 2;
340  std::vector<int> factors;
341  Factorize(N, &factors);
342  int *factor_beg = NULL;
343  if (factors.size() > 0)
344  factor_beg = &(factors[0]);
345  Vector<Real> tmp; // allocated in ComplexFftRecursive.
346  ComplexFftRecursive(v->Data(), 1, N, factor_beg, factor_beg+factors.size(), forward, (tmp_in?tmp_in:&tmp));
347 }
void ComplexFftRecursive(Real *data, int nffts, int N, const int *factor_begin, const int *factor_end, bool forward, Vector< Real > *tmp_vec)
ComplexFftRecursive is a recursive function that computes the complex FFT of size N...
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void Factorize(I m, std::vector< I > *factors)
Definition: kaldi-math.h:322
void ComplexFt ( const VectorBase< Real > &  in,
VectorBase< Real > *  out,
bool  forward 
)

ComplexFt is the same as ComplexFft but it implements the Fourier transform in an inefficient way.

It is mainly included for testing purposes. See comment for ComplexFft to describe the input and outputs and what it does.

Definition at line 29 of file matrix-functions.cc.

References kaldi::ComplexAddProduct(), kaldi::ComplexImExp(), kaldi::ComplexMul(), VectorBase< Real >::Data(), VectorBase< Real >::Dim(), KALDI_ASSERT, and M_2PI.

Referenced by kaldi::UnitTestComplexFft(), kaldi::UnitTestComplexFft2(), kaldi::UnitTestComplexFt(), and kaldi::UnitTestSplitRadixComplexFft().

30  {
31  int exp_sign = (forward ? -1 : 1);
32  KALDI_ASSERT(out != NULL);
33  KALDI_ASSERT(in.Dim() == out->Dim());
34  KALDI_ASSERT(in.Dim() % 2 == 0);
35  int twoN = in.Dim(), N = twoN / 2;
36  const Real *data_in = in.Data();
37  Real *data_out = out->Data();
38 
39  Real exp1N_re, exp1N_im; // forward -> exp(-2pi / N), backward -> exp(2pi / N).
40  Real fraction = exp_sign * M_2PI / static_cast<Real>(N); // forward -> -2pi/N, backward->-2pi/N
41  ComplexImExp(fraction, &exp1N_re, &exp1N_im);
42 
43  Real expm_re = 1.0, expm_im = 0.0; // forward -> exp(-2pi m / N).
44 
45  for (int two_m = 0; two_m < twoN; two_m+=2) { // For each output component.
46  Real expmn_re = 1.0, expmn_im = 0.0; // forward -> exp(-2pi m n / N).
47  Real sum_re = 0.0, sum_im = 0.0; // complex output for index m (the sum expression)
48  for (int two_n = 0; two_n < twoN; two_n+=2) {
49  ComplexAddProduct(data_in[two_n], data_in[two_n+1],
50  expmn_re, expmn_im,
51  &sum_re, &sum_im);
52  ComplexMul(expm_re, expm_im, &expmn_re, &expmn_im);
53  }
54  data_out[two_m] = sum_re;
55  data_out[two_m + 1] = sum_im;
56 
57 
58  if (two_m % 10 == 0) { // occasionally renew "expm" from scratch to avoid
59  // loss of precision.
60  int nextm = 1 + two_m/2;
61  Real fraction_mult = fraction * nextm;
62  ComplexImExp(fraction_mult, &expm_re, &expm_im);
63  } else {
64  ComplexMul(exp1N_re, exp1N_im, &expm_re, &expm_im);
65  }
66  }
67 }
void ComplexAddProduct(const Real &a_re, const Real &a_im, const Real &b_re, const Real &b_im, Real *c_re, Real *c_im)
ComplexMul implements, inline, the complex operation c += (a * b).
void ComplexMul(const Real &a_re, const Real &a_im, Real *b_re, Real *b_im)
ComplexMul implements, inline, the complex multiplication b *= a.
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
#define M_2PI
Definition: kaldi-math.h:52
void ComplexImExp(Real x, Real *a_re, Real *a_im)
ComplexImExp implements a <– exp(i x), inline.
void ComplexImExp ( Real  x,
Real *  a_re,
Real *  a_im 
)
inline

ComplexImExp implements a <– exp(i x), inline.

Definition at line 46 of file matrix-functions-inl.h.

Referenced by kaldi::ComplexFftRecursive(), kaldi::ComplexFt(), SplitRadixRealFft< Real >::Compute(), and kaldi::RealFft().

46  {
47  *a_re = std::cos(x);
48  *a_im = std::sin(x);
49 }
void ComplexMul ( const Real &  a_re,
const Real &  a_im,
Real *  b_re,
Real *  b_im 
)
inline

ComplexMul implements, inline, the complex multiplication b *= a.

Definition at line 31 of file matrix-functions-inl.h.

Referenced by kaldi::ComplexFftRecursive(), kaldi::ComplexFt(), SplitRadixRealFft< Real >::Compute(), kaldi::ElementwiseProductOfFft(), and kaldi::RealFft().

32  {
33  Real tmp_re = (*b_re * a_re) - (*b_im * a_im);
34  *b_im = *b_re * a_im + *b_im * a_re;
35  *b_re = tmp_re;
36 }
void ComputeDctMatrix ( Matrix< Real > *  M)

ComputeDctMatrix computes a matrix corresponding to the DCT, such that M * v equals the DCT of vector v.

M must be square at input. This is the type = III DCT with normalization, corresponding to the following equations, where x is the signal and X is the DCT: X_0 = 1/sqrt(2*N) {n = 0}^{N-1} x_n X_k = 1/sqrt(N) {n = 0}^{N-1} x_n cos( /N (n + 1/2) k ) This matrix's transpose is its own inverse, so transposing this matrix will give the inverse DCT. Caution: the type III DCT is generally known as the "inverse DCT" (with the type II being the actual DCT), so this function is somewhatd mis-named. It was probably done this way for HTK compatibility. We don't change it because it was this way from the start and changing it would affect the feature generation.

Definition at line 592 of file matrix-functions.cc.

References rnnlm::j, KALDI_ASSERT, M_PI, rnnlm::n, MatrixBase< Real >::NumCols(), and MatrixBase< Real >::NumRows().

Referenced by DctComponent::Init(), MfccComputer::MfccComputer(), and kaldi::UnitTestDct().

592  {
593  //KALDI_ASSERT(M->NumRows() == M->NumCols());
594  MatrixIndexT K = M->NumRows();
595  MatrixIndexT N = M->NumCols();
596 
597  KALDI_ASSERT(K > 0);
598  KALDI_ASSERT(N > 0);
599  Real normalizer = std::sqrt(1.0 / static_cast<Real>(N)); // normalizer for
600  // X_0.
601  for (MatrixIndexT j = 0; j < N; j++) (*M)(0, j) = normalizer;
602  normalizer = std::sqrt(2.0 / static_cast<Real>(N)); // normalizer for other
603  // elements.
604  for (MatrixIndexT k = 1; k < K; k++)
605  for (MatrixIndexT n = 0; n < N; n++)
606  (*M)(k, n) = normalizer
607  * std::cos( static_cast<double>(M_PI)/N * (n + 0.5) * k );
608 }
#define M_PI
Definition: kaldi-math.h:44
int32 MatrixIndexT
Definition: matrix-common.h:96
struct rnnlm::@11::@12 n
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void ComputePca ( const MatrixBase< Real > &  X,
MatrixBase< Real > *  U,
MatrixBase< Real > *  A,
bool  print_eigs = false,
bool  exact = true 
)

ComputePCA does a PCA computation, using either outer products or inner products, whichever is more efficient.

Let D be the dimension of the data points, N be the number of data points, and G be the PCA dimension we want to retain. We assume G <= N and G <= D.

Parameters
X[in] An N x D matrix. Each row of X is a point x_i.
U[out] A G x D matrix. Each row of U is a basis element u_i.
A[out] An N x D matrix, or NULL. Each row of A is a set of coefficients in the basis for a point x_i, so A(i, g) is the coefficient of u_i in x_i.
print_eigs[in] If true, prints out diagnostic information about the eigenvalues.
exact[in] If true, does the exact computation; if false, does a much faster (but almost exact) computation based on the Lanczos method.

Definition at line 823 of file matrix-functions.cc.

References SpMatrix< Real >::AddMat2(), MatrixBase< Real >::AddMatMat(), MatrixBase< Real >::DestructiveSvd(), SpMatrix< Real >::Eig(), KALDI_ASSERT, KALDI_LOG, KALDI_WARN, kaldi::kNoTrans, kaldi::kTrans, rnnlm::n, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), MatrixBase< Real >::OrthogonalizeRows(), Vector< Real >::Resize(), Matrix< Real >::Resize(), MatrixBase< Real >::Row(), kaldi::SortSvd(), SpMatrix< Real >::TopEigs(), and Matrix< Real >::Transpose().

Referenced by kaldi::UnitTestPca(), and kaldi::UnitTestPca2().

827  {
828  // Note that some of these matrices may be transposed w.r.t. the
829  // way it's most natural to describe them in math... it's the rows
830  // of X and U that correspond to the (data-points, basis elements).
831  MatrixIndexT N = X.NumRows(), D = X.NumCols();
832  // N = #points, D = feature dim.
833  KALDI_ASSERT(U != NULL && U->NumCols() == D);
834  MatrixIndexT G = U->NumRows(); // # of retained basis elements.
835  KALDI_ASSERT(A == NULL || (A->NumRows() == N && A->NumCols() == G));
836  KALDI_ASSERT(G <= N && G <= D);
837  if (D < N) { // Do conventional PCA.
838  SpMatrix<Real> Msp(D); // Matrix of outer products.
839  Msp.AddMat2(1.0, X, kTrans, 0.0); // M <-- X^T X
840  Matrix<Real> Utmp;
841  Vector<Real> l;
842  if (exact) {
843  Utmp.Resize(D, D);
844  l.Resize(D);
845  //Matrix<Real> M(Msp);
846  //M.DestructiveSvd(&l, &Utmp, NULL);
847  Msp.Eig(&l, &Utmp);
848  } else {
849  Utmp.Resize(D, G);
850  l.Resize(G);
851  Msp.TopEigs(&l, &Utmp);
852  }
853  SortSvd(&l, &Utmp);
854 
855  for (MatrixIndexT g = 0; g < G; g++)
856  U->Row(g).CopyColFromMat(Utmp, g);
857  if (print_eigs)
858  KALDI_LOG << (exact ? "" : "Retained ")
859  << "PCA eigenvalues are " << l;
860  if (A != NULL)
861  A->AddMatMat(1.0, X, kNoTrans, *U, kTrans, 0.0);
862  } else { // Do inner-product PCA.
863  SpMatrix<Real> Nsp(N); // Matrix of inner products.
864  Nsp.AddMat2(1.0, X, kNoTrans, 0.0); // M <-- X X^T
865 
866  Matrix<Real> Vtmp;
867  Vector<Real> l;
868  if (exact) {
869  Vtmp.Resize(N, N);
870  l.Resize(N);
871  Matrix<Real> Nmat(Nsp);
872  Nmat.DestructiveSvd(&l, &Vtmp, NULL);
873  } else {
874  Vtmp.Resize(N, G);
875  l.Resize(G);
876  Nsp.TopEigs(&l, &Vtmp);
877  }
878 
879  MatrixIndexT num_zeroed = 0;
880  for (MatrixIndexT g = 0; g < G; g++) {
881  if (l(g) < 0.0) {
882  KALDI_WARN << "In PCA, setting element " << l(g) << " to zero.";
883  l(g) = 0.0;
884  num_zeroed++;
885  }
886  }
887  SortSvd(&l, &Vtmp); // Make sure zero elements are last, this
888  // is necessary for Orthogonalize() to work properly later.
889 
890  Vtmp.Transpose(); // So eigenvalues are the rows.
891 
892  for (MatrixIndexT g = 0; g < G; g++) {
893  Real sqrtlg = sqrt(l(g));
894  if (l(g) != 0.0) {
895  U->Row(g).AddMatVec(1.0 / sqrtlg, X, kTrans, Vtmp.Row(g), 0.0);
896  } else {
897  U->Row(g).SetZero();
898  (*U)(g, g) = 1.0; // arbitrary direction. Will later orthogonalize.
899  }
900  if (A != NULL)
901  for (MatrixIndexT n = 0; n < N; n++)
902  (*A)(n, g) = sqrtlg * Vtmp(g, n);
903  }
904  // Now orthogonalize. This is mainly useful in
905  // case there were zero eigenvalues, but we do it
906  // for all of them.
907  U->OrthogonalizeRows();
908  if (print_eigs)
909  KALDI_LOG << "(inner-product) PCA eigenvalues are " << l;
910  }
911 }
int32 MatrixIndexT
Definition: matrix-common.h:96
struct rnnlm::@11::@12 n
#define KALDI_WARN
Definition: kaldi-error.h:130
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
#define KALDI_LOG
Definition: kaldi-error.h:133
void SortSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt, bool sort_on_absolute_value)
Function to ensure that SVD is sorted.
void CreateEigenvalueMatrix ( const VectorBase< Real > &  real,
const VectorBase< Real > &  imag,
MatrixBase< Real > *  D 
)

Creates the eigenvalue matrix D that is part of the decomposition used Matrix::Eig.

D will be block-diagonal with blocks of size 1 (for real eigenvalues) or 2x2 for complex pairs. If a complex pair is lambda +- i*mu, D will have a corresponding 2x2 block [lambda, mu; -mu, lambda]. This function will throw if any complex eigenvalues are not in complex conjugate pairs (or the members of such pairs are not consecutively numbered).

if (im(j) < 0.0) KALDI_WARN << "Negative first im part of pair"; // TEMP

Definition at line 2411 of file kaldi-matrix.cc.

References kaldi::ApproxEqual(), VectorBase< Real >::Dim(), rnnlm::j, KALDI_ASSERT, rnnlm::n, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), and MatrixBase< Real >::SetZero().

Referenced by MatrixBase< Real >::Power(), and kaldi::UnitTestEig().

2412  {
2413  MatrixIndexT n = re.Dim();
2414  KALDI_ASSERT(im.Dim() == n && D->NumRows() == n && D->NumCols() == n);
2415 
2416  MatrixIndexT j = 0;
2417  D->SetZero();
2418  while (j < n) {
2419  if (im(j) == 0) { // Real eigenvalue
2420  (*D)(j, j) = re(j);
2421  j++;
2422  } else { // First of a complex pair
2423  KALDI_ASSERT(j+1 < n && ApproxEqual(im(j+1), -im(j))
2424  && ApproxEqual(re(j+1), re(j)));
2426  Real lambda = re(j), mu = im(j);
2427  // create 2x2 block [lambda, mu; -mu, lambda]
2428  (*D)(j, j) = lambda;
2429  (*D)(j, j+1) = mu;
2430  (*D)(j+1, j) = -mu;
2431  (*D)(j+1, j+1) = lambda;
2432  j += 2;
2433  }
2434  }
2435 }
int32 MatrixIndexT
Definition: matrix-common.h:96
struct rnnlm::@11::@12 n
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
static bool ApproxEqual(float a, float b, float relative_tolerance=0.001)
return abs(a - b) <= relative_tolerance * (abs(a)+abs(b)).
Definition: kaldi-math.h:262
void RealFft ( VectorBase< Real > *  v,
bool  forward 
)

RealFft is a fourier transform of real inputs.

Internally it uses ComplexFft. The input dimension N must be even. If forward == true, it transforms from a sequence of N real points to its complex fourier transform; otherwise it goes in the reverse direction. If you call it in the forward and then reverse direction and multiply by 1.0/N, you will get back the original data. The interpretation of the complex-FFT data is as follows: the array is a sequence of complex numbers C_n of length N/2 with (real, im) format, i.e. [real0, real_{N/2}, real1, im1, real2, im2, real3, im3, ...]. See also SplitRadixRealFft, declared in srfft.h, which is more efficient but only works if the length of the input is a power of 2.

Definition at line 391 of file matrix-functions.cc.

References kaldi::ComplexAddProduct(), kaldi::ComplexFft(), kaldi::ComplexImExp(), kaldi::ComplexMul(), VectorBase< Real >::Data(), VectorBase< Real >::Dim(), KALDI_ASSERT, M_2PI, and VectorBase< Real >::Scale().

Referenced by SpectrogramComputer::Compute(), MfccComputer::Compute(), FbankComputer::Compute(), PlpComputer::Compute(), kaldi::UnitTestRealFft(), and kaldi::UnitTestRealFftSpeed().

391  {
392  KALDI_ASSERT(v != NULL);
393  MatrixIndexT N = v->Dim(), N2 = N/2;
394  KALDI_ASSERT(N%2 == 0);
395  if (N == 0) return;
396 
397  if (forward) ComplexFft(v, true);
398 
399  Real *data = v->Data();
400  Real rootN_re, rootN_im; // exp(-2pi/N), forward; exp(2pi/N), backward
401  int forward_sign = forward ? -1 : 1;
402  ComplexImExp(static_cast<Real>(M_2PI/N *forward_sign), &rootN_re, &rootN_im);
403  Real kN_re = -forward_sign, kN_im = 0.0; // exp(-2pik/N), forward; exp(-2pik/N), backward
404  // kN starts out as 1.0 for forward algorithm but -1.0 for backward.
405  for (MatrixIndexT k = 1; 2*k <= N2; k++) {
406  ComplexMul(rootN_re, rootN_im, &kN_re, &kN_im);
407 
408  Real Ck_re, Ck_im, Dk_re, Dk_im;
409  // C_k = 1/2 (B_k + B_{N/2 - k}^*) :
410  Ck_re = 0.5 * (data[2*k] + data[N - 2*k]);
411  Ck_im = 0.5 * (data[2*k + 1] - data[N - 2*k + 1]);
412  // re(D_k)= 1/2 (im(B_k) + im(B_{N/2-k})):
413  Dk_re = 0.5 * (data[2*k + 1] + data[N - 2*k + 1]);
414  // im(D_k) = -1/2 (re(B_k) - re(B_{N/2-k}))
415  Dk_im =-0.5 * (data[2*k] - data[N - 2*k]);
416  // A_k = C_k + 1^(k/N) D_k:
417  data[2*k] = Ck_re; // A_k <-- C_k
418  data[2*k+1] = Ck_im;
419  // now A_k += D_k 1^(k/N)
420  ComplexAddProduct(Dk_re, Dk_im, kN_re, kN_im, &(data[2*k]), &(data[2*k+1]));
421 
422  MatrixIndexT kdash = N2 - k;
423  if (kdash != k) {
424  // Next we handle the index k' = N/2 - k. This is necessary
425  // to do now, to avoid invalidating data that we will later need.
426  // The quantities C_{k'} and D_{k'} are just the conjugates of C_k
427  // and D_k, so the equations are simple modifications of the above,
428  // replacing Ck_im and Dk_im with their negatives.
429  data[2*kdash] = Ck_re; // A_k' <-- C_k'
430  data[2*kdash+1] = -Ck_im;
431  // now A_k' += D_k' 1^(k'/N)
432  // We use 1^(k'/N) = 1^((N/2 - k) / N) = 1^(1/2) 1^(-k/N) = -1 * (1^(k/N))^*
433  // so it's the same as 1^(k/N) but with the real part negated.
434  ComplexAddProduct(Dk_re, -Dk_im, -kN_re, kN_im, &(data[2*kdash]), &(data[2*kdash+1]));
435  }
436  }
437 
438  { // Now handle k = 0.
439  // In simple terms: after the complex fft, data[0] becomes the sum of real
440  // parts input[0], input[2]... and data[1] becomes the sum of imaginary
441  // pats input[1], input[3]...
442  // "zeroth" [A_0] is just the sum of input[0]+input[1]+input[2]..
443  // and "n2th" [A_{N/2}] is input[0]-input[1]+input[2]... .
444  Real zeroth = data[0] + data[1],
445  n2th = data[0] - data[1];
446  data[0] = zeroth;
447  data[1] = n2th;
448  if (!forward) {
449  data[0] /= 2;
450  data[1] /= 2;
451  }
452  }
453 
454  if (!forward) {
455  ComplexFft(v, false);
456  v->Scale(2.0); // This is so we get a factor of N increase, rather than N/2 which we would
457  // otherwise get from [ComplexFft, forward] + [ComplexFft, backward] in dimension N/2.
458  // It's for consistency with our normal FFT convensions.
459  }
460 }
template void ComplexFft(VectorBase< double > *v, bool forward, Vector< double > *tmp_in)
int32 MatrixIndexT
Definition: matrix-common.h:96
void ComplexAddProduct(const Real &a_re, const Real &a_im, const Real &b_re, const Real &b_im, Real *c_re, Real *c_im)
ComplexMul implements, inline, the complex operation c += (a * b).
void ComplexMul(const Real &a_re, const Real &a_im, Real *b_re, Real *b_im)
ComplexMul implements, inline, the complex multiplication b *= a.
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
#define M_2PI
Definition: kaldi-math.h:52
void ComplexImExp(Real x, Real *a_re, Real *a_im)
ComplexImExp implements a <– exp(i x), inline.
void RealFftInefficient ( VectorBase< Real > *  v,
bool  forward 
)

Inefficient version of Fourier transform, for testing purposes.

RealFt has the same input and output format as RealFft above, but it is an inefficient implementation included for testing purposes.

Definition at line 350 of file matrix-functions.cc.

References kaldi::ComplexFft(), VectorBase< Real >::CopyFromVec(), VectorBase< Real >::Dim(), rnnlm::i, KALDI_ASSERT, and VectorBase< Real >::Range().

Referenced by kaldi::UnitTestRealFft(), and kaldi::UnitTestSplitRadixRealFft().

350  {
351  KALDI_ASSERT(v != NULL);
352  MatrixIndexT N = v->Dim();
353  KALDI_ASSERT(N%2 == 0);
354  if (N == 0) return;
355  Vector<Real> vtmp(N*2); // store as complex.
356  if (forward) {
357  for (MatrixIndexT i = 0; i < N; i++) vtmp(i*2) = (*v)(i);
358  ComplexFft(&vtmp, forward); // this is already tested so we can use this.
359  v->CopyFromVec( vtmp.Range(0, N) );
360  (*v)(1) = vtmp(N); // Copy the N/2'th fourier component, which is real,
361  // to the imaginary part of the 1st complex output.
362  } else {
363  // reverse the transformation above to get the complex spectrum.
364  vtmp(0) = (*v)(0); // copy F_0 which is real
365  vtmp(N) = (*v)(1); // copy F_{N/2} which is real
366  for (MatrixIndexT i = 1; i < N/2; i++) {
367  // Copy i'th to i'th fourier component
368  vtmp(2*i) = (*v)(2*i);
369  vtmp(2*i+1) = (*v)(2*i+1);
370  // Copy i'th to N-i'th, conjugated.
371  vtmp(2*(N-i)) = (*v)(2*i);
372  vtmp(2*(N-i)+1) = -(*v)(2*i+1);
373  }
374  ComplexFft(&vtmp, forward); // actually backward since forward == false
375  // Copy back real part. Complex part should be zero.
376  for (MatrixIndexT i = 0; i < N; i++)
377  (*v)(i) = vtmp(i*2);
378  }
379 }
template void ComplexFft(VectorBase< double > *v, bool forward, Vector< double > *tmp_in)
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
Real SolveDoubleQuadraticMatrixProblem ( const MatrixBase< Real > &  G,
const SpMatrix< Real > &  P1,
const SpMatrix< Real > &  P2,
const SpMatrix< Real > &  Q1,
const SpMatrix< Real > &  Q2,
const SolverOptions &  opts,
MatrixBase< Real > *  M 
)

Maximizes the auxiliary function :

\[ Q(M) = tr(M^T G) -0.5 tr(P_1 M Q_1 M^T) -0.5 tr(P_2 M Q_2 M^T). \]

Encountered in matrix update with a prior.

We also apply a limit on the condition but it should be less frequently necessary, and can be set larger.

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

References SpMatrix< Real >::AddMat2Sp(), MatrixBase< Real >::AddMatMat(), VectorBase< Real >::AddMatVec(), SpMatrix< Real >::AddSp(), TpMatrix< Real >::Cholesky(), rnnlm::d, TpMatrix< Real >::Invert(), MatrixBase< Real >::Invert(), SpMatrix< Real >::IsDiagonal(), SpMatrix< Real >::IsUnit(), KALDI_ASSERT, KALDI_ERR, KALDI_WARN, kaldi::kNoTrans, kaldi::kTrans, rnnlm::n, SolverOptions::name, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), PackedMatrix< Real >::NumRows(), MatrixBase< Real >::Row(), kaldi::SolveQuadraticProblem(), SpMatrix< Real >::SymPosSemiDefEig(), kaldi::VecSpVec(), and kaldi::VecVec().

Referenced by MleAmSgmm2Updater::MapUpdateM(), and kaldi::UnitTestSolve().

857  {
858  KALDI_ASSERT(Q1.NumRows() == M->NumCols() && P1.NumRows() == M->NumRows() &&
859  G.NumRows() == M->NumRows() && G.NumCols() == M->NumCols() &&
860  M->NumCols() != 0 && Q2.NumRows() == M->NumCols() &&
861  P2.NumRows() == M->NumRows());
862  MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
863  // The following check should not fail as we stipulate P1, P2 and one of Q1
864  // or Q2 must be +ve def and other Q1 or Q2 must be +ve semidef.
865  TpMatrix<Real> LInv(rows);
866  LInv.Cholesky(P1);
867  LInv.Invert(); // Will throw exception if fails.
868  SpMatrix<Real> S(rows);
869  Matrix<Real> LInvFull(LInv);
870  S.AddMat2Sp(1.0, LInvFull, kNoTrans, P2, 0.0); // S := L^{-1} P_2 L^{-T}
871  Matrix<Real> U(rows, rows);
872  Vector<Real> d(rows);
873  S.SymPosSemiDefEig(&d, &U);
874  Matrix<Real> T(rows, rows);
875  T.AddMatMat(1.0, U, kTrans, LInvFull, kNoTrans, 0.0); // T := U^T * L^{-1}
876 
877 #ifdef KALDI_PARANOID // checking mainly for errors in the code or math.
878  {
879  SpMatrix<Real> P1Trans(rows);
880  P1Trans.AddMat2Sp(1.0, T, kNoTrans, P1, 0.0);
881  KALDI_ASSERT(P1Trans.IsUnit(0.01));
882  }
883  {
884  SpMatrix<Real> P2Trans(rows);
885  P2Trans.AddMat2Sp(1.0, T, kNoTrans, P2, 0.0);
886  KALDI_ASSERT(P2Trans.IsDiagonal(0.01));
887  }
888 #endif
889 
890  Matrix<Real> TInv(T);
891  TInv.Invert();
892  Matrix<Real> Gdash(rows, cols);
893  Gdash.AddMatMat(1.0, T, kNoTrans, G, kNoTrans, 0.0); // G' = T G
894  Matrix<Real> MdashOld(rows, cols);
895  MdashOld.AddMatMat(1.0, TInv, kTrans, *M, kNoTrans, 0.0); // M' = T^{-T} M
896  Matrix<Real> MdashNew(MdashOld);
897  Real objf_impr = 0.0;
898  for (MatrixIndexT n = 0; n < rows; n++) {
899  SpMatrix<Real> Qsum(Q1);
900  Qsum.AddSp(d(n), Q2);
901  SubVector<Real> mdash_n = MdashNew.Row(n);
902  SubVector<Real> gdash_n = Gdash.Row(n);
903 
904  Matrix<Real> QsumInv(Qsum);
905  try {
906  QsumInv.Invert();
907  Real old_objf = VecVec(mdash_n, gdash_n)
908  - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
909  mdash_n.AddMatVec(1.0, QsumInv, kNoTrans, gdash_n, 0.0); // m'_n := g'_n * (Q_1 + d_n Q_2)^{-1}
910  Real new_objf = VecVec(mdash_n, gdash_n)
911  - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
912  if (new_objf < old_objf) {
913  if (new_objf < old_objf - 1.0e-05) {
914  KALDI_WARN << "In double quadratic matrix problem: objective "
915  "function decreasing during optimization of " << opts.name
916  << ", " << old_objf << "->" << new_objf << ", change is "
917  << (new_objf - old_objf);
918  KALDI_ERR << "Auxiliary function decreasing."; // Will be caught.
919  } else { // Reset to old value, didn't improve (very close to optimum).
920  MdashNew.Row(n).CopyFromVec(MdashOld.Row(n));
921  }
922  }
923  objf_impr += new_objf - old_objf;
924  }
925  catch (...) {
926  KALDI_WARN << "Matrix inversion or optimization failed during double "
927  "quadratic problem, solving for" << opts.name
928  << ": trying more stable approach.";
929  objf_impr += SolveQuadraticProblem(Qsum, gdash_n, opts, &mdash_n);
930  }
931  }
932  M->AddMatMat(1.0, T, kTrans, MdashNew, kNoTrans, 0.0); // M := T^T M'.
933  return objf_impr;
934 }
double SolveQuadraticProblem(const SpMatrix< double > &H, const VectorBase< double > &g, const SolverOptions &opts, VectorBase< double > *x)
Definition: sp-matrix.cc:659
template double VecSpVec(const VectorBase< double > &v1, const SpMatrix< double > &M, const VectorBase< double > &v2)
int32 MatrixIndexT
Definition: matrix-common.h:96
struct rnnlm::@11::@12 n
#define KALDI_ERR
Definition: kaldi-error.h:127
#define KALDI_WARN
Definition: kaldi-error.h:130
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
Real VecVec(const VectorBase< Real > &a, const VectorBase< Real > &b)
Returns dot product between v1 and v2.
Definition: kaldi-vector.cc:36
Real SolveQuadraticMatrixProblem ( const SpMatrix< Real > &  Q,
const MatrixBase< Real > &  Y,
const SpMatrix< Real > &  P,
const SolverOptions &  opts,
MatrixBase< Real > *  M 
)

Maximizes the auxiliary function :

\[ Q(x) = tr(M^T P Y) - 0.5 tr(P M Q M^T) \]

Like a numerically stable version of $ M := Y Q^{-1} $.

Assumes Q and P positive semidefinite, and matrix dimensions match enough to make expressions meaningful. This is mostly as described in the SGMM paper "the subspace Gaussian mixture model &ndash; a structured model for speech recognition", but the diagonal_precondition option is newly added, to handle problems where different dimensions have very different scaling (we recommend to use the option but it's set false for back compatibility).

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

References MatrixBase< Real >::AddMat(), SpMatrix< Real >::AddMat2Sp(), MatrixBase< Real >::AddMatMat(), SpMatrix< Real >::AddVec2Sp(), VectorBase< Real >::ApplyFloor(), VectorBase< Real >::ApplyPow(), SolverOptions::Check(), VectorBase< Real >::CopyDiagFromSp(), MatrixBase< Real >::CopyFromMat(), SolverOptions::diagonal_precondition, SolverOptions::eps, rnnlm::i, VectorBase< Real >::InvertElements(), SpMatrix< Real >::IsZero(), SolverOptions::K, KALDI_ASSERT, KALDI_LOG, KALDI_WARN, kaldi::kNoTrans, kaldi::kTrans, VectorBase< Real >::Max(), MatrixBase< Real >::MulColsVec(), SolverOptions::name, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), PackedMatrix< Real >::NumRows(), SolverOptions::optimize_delta, SolverOptions::print_debug_output, SpMatrix< Real >::SymPosSemiDefEig(), kaldi::TraceMatMat(), and kaldi::TraceSpSp().

Referenced by kaldi::UnitTestSolve(), EbwAmSgmm2Updater::UpdateM(), MleAmSgmm2Updater::UpdateM(), EbwAmSgmm2Updater::UpdateN(), MleAmSgmm2Updater::UpdateN(), and IvectorExtractorStats::UpdateProjection().

757  {
758  KALDI_ASSERT(Q.NumRows() == M->NumCols() &&
759  SigmaInv.NumRows() == M->NumRows() && Y.NumRows() == M->NumRows()
760  && Y.NumCols() == M->NumCols() && M->NumCols() != 0);
761  opts.Check();
762  MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
763  if (Q.IsZero(0.0)) {
764  KALDI_WARN << "Zero quadratic term in quadratic matrix problem for "
765  << opts.name << ": leaving it unchanged.";
766  return 0.0;
767  }
768 
769  if (opts.diagonal_precondition) {
770  // We can re-cast the problem with a diagonal preconditioner in the space
771  // of Q (columns of M). Helps to improve the condition of Q.
772  Vector<Real> Q_diag(cols);
773  Q_diag.CopyDiagFromSp(Q);
774  Q_diag.ApplyFloor(std::numeric_limits<Real>::min() * 1.0E+3);
775  Vector<Real> Q_diag_sqrt(Q_diag);
776  Q_diag_sqrt.ApplyPow(0.5);
777  Vector<Real> Q_diag_inv_sqrt(Q_diag_sqrt);
778  Q_diag_inv_sqrt.InvertElements();
779  Matrix<Real> M_scaled(*M);
780  M_scaled.MulColsVec(Q_diag_sqrt);
781  Matrix<Real> Y_scaled(Y);
782  Y_scaled.MulColsVec(Q_diag_inv_sqrt);
783  SpMatrix<Real> Q_scaled(cols);
784  Q_scaled.AddVec2Sp(1.0, Q_diag_inv_sqrt, Q, 0.0);
785  Real ans;
786  SolverOptions new_opts(opts);
787  new_opts.diagonal_precondition = false;
788  ans = SolveQuadraticMatrixProblem(Q_scaled, Y_scaled, SigmaInv,
789  new_opts, &M_scaled);
790  M->CopyFromMat(M_scaled);
791  M->MulColsVec(Q_diag_inv_sqrt);
792  return ans;
793  }
794 
795  Matrix<Real> Ybar(Y);
796  if (opts.optimize_delta) {
797  Matrix<Real> Qfull(Q);
798  Ybar.AddMatMat(-1.0, *M, kNoTrans, Qfull, kNoTrans, 1.0);
799  } // Ybar = Y - M Q.
800  Matrix<Real> U(cols, cols);
801  Vector<Real> l(cols);
802  Q.SymPosSemiDefEig(&l, &U); // does svd Q = U L V^T and checks that Q == U L U^T to within a tolerance.
803  // floor l.
804  Real f = std::max<Real>(static_cast<Real>(opts.eps), l.Max() / opts.K);
805  MatrixIndexT nfloored = 0;
806  for (MatrixIndexT i = 0; i < cols; i++) { // floor l.
807  if (l(i) < f) { nfloored++; l(i) = f; }
808  }
809  if (nfloored != 0 && opts.print_debug_output)
810  KALDI_LOG << "Solving matrix problem for " << opts.name
811  << ": floored " << nfloored << " eigenvalues. ";
812  Matrix<Real> tmpDelta(rows, cols);
813  tmpDelta.AddMatMat(1.0, Ybar, kNoTrans, U, kNoTrans, 0.0); // tmpDelta = Ybar * U.
814  l.InvertElements(); KALDI_ASSERT(1.0/l.Max() != 0); // check not infinite. eps should take care of this.
815  tmpDelta.MulColsVec(l); // tmpDelta = Ybar * U * \tilde{L}^{-1}
816 
817  Matrix<Real> Delta(rows, cols);
818  Delta.AddMatMat(1.0, tmpDelta, kNoTrans, U, kTrans, 0.0); // Delta = Ybar * U * \tilde{L}^{-1} * U^T
819 
820  Real auxf_before, auxf_after;
821  SpMatrix<Real> MQM(rows);
822  Matrix<Real> &SigmaInvY(tmpDelta);
823  { Matrix<Real> SigmaInvFull(SigmaInv); SigmaInvY.AddMatMat(1.0, SigmaInvFull, kNoTrans, Y, kNoTrans, 0.0); }
824  { // get auxf_before. Q(x) = tr(M^T SigmaInv Y) - 0.5 tr(SigmaInv M Q M^T).
825  MQM.AddMat2Sp(1.0, *M, kNoTrans, Q, 0.0);
826  auxf_before = TraceMatMat(*M, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
827  }
828 
829  Matrix<Real> Mhat(Delta);
830  if (opts.optimize_delta) Mhat.AddMat(1.0, *M); // Mhat = Delta + M.
831 
832  { // get auxf_after.
833  MQM.AddMat2Sp(1.0, Mhat, kNoTrans, Q, 0.0);
834  auxf_after = TraceMatMat(Mhat, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
835  }
836 
837  if (auxf_after < auxf_before) {
838  if (auxf_after < auxf_before - 1.0e-10)
839  KALDI_WARN << "Optimizing matrix auxiliary function for "
840  << opts.name << ", auxf decreased "
841  << auxf_before << " to " << auxf_after << ", change is "
842  << (auxf_after-auxf_before);
843  return 0.0;
844  } else {
845  M->CopyFromMat(Mhat);
846  return auxf_after - auxf_before;
847  }
848 }
int32 MatrixIndexT
Definition: matrix-common.h:96
template double SolveQuadraticMatrixProblem(const SpMatrix< double > &Q, const MatrixBase< double > &Y, const SpMatrix< double > &SigmaInv, const SolverOptions &opts, MatrixBase< double > *M)
#define KALDI_WARN
Definition: kaldi-error.h:130
Real TraceMatMat(const MatrixBase< Real > &A, const MatrixBase< Real > &B, MatrixTransposeType trans)
We need to declare this here as it will be a friend function.
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
Real TraceSpSp(const SpMatrix< Real > &A, const SpMatrix< OtherReal > &B)
Returns tr(A B).
Definition: sp-matrix.cc:391
#define KALDI_LOG
Definition: kaldi-error.h:133
Real kaldi::SolveQuadraticProblem ( const SpMatrix< Real > &  H,
const VectorBase< Real > &  g,
const SolverOptions &  opts,
VectorBase< Real > *  x 
)

Maximizes the auxiliary function

\[ Q(x) = x.g - 0.5 x^T H x \]

using a numerically stable method.

Like a numerically stable version of $ x := Q^{-1} g. $ Assumes H positive semidefinite. Returns the objective-function change.

void SortSvd ( VectorBase< Real > *  s,
MatrixBase< Real > *  U,
MatrixBase< Real > *  Vt = NULL,
bool  sort_on_absolute_value = true 
)

Function to ensure that SVD is sorted.

This function is made as generic as possible, to be applicable to other types of problems. s->Dim() should be the same as U->NumCols(), and we sort s from greatest to least absolute value (if sort_on_absolute_value == true) or greatest to least value otherwise, moving the columns of U, if it exists, and the rows of Vt, if it exists, around in the same way. Note: the "absolute value" part won't matter if this is an actual SVD, since singular values are non-negative.

Makes sure the Svd is sorted (from greatest to least absolute value).

Definition at line 2368 of file kaldi-matrix.cc.

References rnnlm::d, VectorBase< Real >::Dim(), KALDI_ASSERT, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), and MatrixBase< Real >::Row().

Referenced by kaldi::ComputeFeatureNormalizingTransform(), kaldi::ComputeLdaTransform(), kaldi::ComputePca(), LdaEstimate::Estimate(), BasisFmllrEstimate::EstimateFmllrBasis(), FeatureTransformEstimate::EstimateInternal(), kaldi::EstimateSgmm2FmllrSubspace(), IvectorExtractorStats::GetOrthogonalIvectorTransform(), PldaEstimator::GetOutput(), AffineComponent::LimitRank(), main(), LimitRankClass::operator()(), OnlinePreconditionerSimple::PreconditionDirectionsCpu(), OnlineNaturalGradientSimple::PreconditionDirectionsCpu(), OnlinePreconditioner::PreconditionDirectionsInternal(), OnlineNaturalGradient::PreconditionDirectionsInternal(), SpMatrix< Real >::TopEigs(), kaldi::UnitTestPca2(), kaldi::UnitTestSvdNodestroy(), and PldaUnsupervisedAdaptor::UpdatePlda().

2369  {
2371  MatrixIndexT num_singval = s->Dim();
2372  KALDI_ASSERT(U == NULL || U->NumCols() == num_singval);
2373  KALDI_ASSERT(Vt == NULL || Vt->NumRows() == num_singval);
2374 
2375  std::vector<std::pair<Real, MatrixIndexT> > vec(num_singval);
2376  // negative because we want revese order.
2377  for (MatrixIndexT d = 0; d < num_singval; d++) {
2378  Real val = (*s)(d),
2379  sort_val = -(sort_on_absolute_value ? std::abs(val) : val);
2380  vec[d] = std::pair<Real, MatrixIndexT>(sort_val, d);
2381  }
2382  std::sort(vec.begin(), vec.end());
2383  Vector<Real> s_copy(*s);
2384  for (MatrixIndexT d = 0; d < num_singval; d++)
2385  (*s)(d) = s_copy(vec[d].second);
2386  if (U != NULL) {
2387  Matrix<Real> Utmp(*U);
2388  MatrixIndexT dim = Utmp.NumRows();
2389  for (MatrixIndexT d = 0; d < num_singval; d++) {
2390  MatrixIndexT oldidx = vec[d].second;
2391  for (MatrixIndexT e = 0; e < dim; e++)
2392  (*U)(e, d) = Utmp(e, oldidx);
2393  }
2394  }
2395  if (Vt != NULL) {
2396  Matrix<Real> Vttmp(*Vt);
2397  for (MatrixIndexT d = 0; d < num_singval; d++)
2398  (*Vt).Row(d).CopyFromVec(Vttmp.Row(vec[d].second));
2399  }
2400 }
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169