Miscellaneous matrix/vector library functions and classes
Collaboration diagram for Miscellaneous matrix/vector library functions and classes:

## Classes

struct  SolverOptions
This class describes the options for maximizing various quadratic objective functions. More...

## 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

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 :

Like a numerically stable version of . 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 :

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

## Function Documentation

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

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

730  {
731  KALDI_ASSERT(a.Dim() == plus->NumRows() && b.Dim() == plus->NumCols()
732  && a.Dim() == minus->NumRows() && b.Dim() == minus->NumCols());
733  int32 nrows = a.Dim(), ncols = b.Dim(), pskip = plus->Stride() - ncols,
734  mskip = minus->Stride() - ncols;
735  const Real *adata = a.Data(), *bdata = b.Data();
736  Real *plusdata = plus->Data(), *minusdata = minus->Data();
737
738  for (int32 i = 0; i < nrows; i++) {
739  const Real *btmp = bdata;
740  Real multiple = alpha * *adata;
741  if (multiple > 0.0) {
742  for (int32 j = 0; j < ncols; j++, plusdata++, minusdata++, btmp++) {
743  if (*btmp > 0.0) *plusdata += multiple * *btmp;
744  else *minusdata -= multiple * *btmp;
745  }
746  } else {
747  for (int32 j = 0; j < ncols; j++, plusdata++, minusdata++, btmp++) {
748  if (*btmp < 0.0) *plusdata += multiple * *btmp;
749  else *minusdata -= multiple * *btmp;
750  }
751  }
752  plusdata += pskip;
753  minusdata += mskip;
755  }
756 }
kaldi::int32 int32
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185

## ◆ AssertSameDim()

 void kaldi::AssertSameDim ( const MatrixBase< Real1 > & mat1, const MatrixBase< Real2 > & mat2 )
inline

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

161  {
162  KALDI_ASSERT(mat1.NumRows() == mat2.NumRows()
163  && mat1.NumCols() == mat2.NumCols());
164 }
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185

## ◆ AttemptComplexPower()

 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 2659 of file kaldi-matrix.cc.

References kaldi::AttemptComplexPower().

Referenced by kaldi::TraceMat(), and kaldi::UnitTestComplexPower().

2659  {
2660  // Used in Matrix<Real>::Power().
2661  // Attempts to take the complex value x to the power "power",
2662  // assuming that power is fractional (i.e. we don't treat integers as a
2663  // special case). Returns false if this is not possible, either
2664  // because x is negative and real (hence there is no obvious answer
2665  // that is "closest to 1", and anyway this case does not make sense
2666  // in the Matrix<Real>::Power() routine);
2667  // or because power is negative, and x is zero.
2668
2669  // First solve for r and theta in
2670  // x_re = r*cos(theta), x_im = r*sin(theta)
2671  if (*x_re < 0.0 && *x_im == 0.0) return false; // can't do
2672  // it for negative real values.
2673  Real r = std::sqrt((*x_re * *x_re) + (*x_im * *x_im)); // r == radius.
2674  if (power < 0.0 && r == 0.0) return false;
2675  Real theta = std::atan2(*x_im, *x_re);
2676  // Take the power.
2677  r = std::pow(r, power);
2678  theta *= power;
2679  *x_re = r * std::cos(theta);
2680  *x_im = r * std::sin(theta);
2681  return true;
2682 }

 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.

40  {
41  *c_re += b_re*a_re - b_im*a_im;
42  *c_im += b_re*a_im + b_im*a_re;
43 }

## ◆ ComplexFft()

 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.

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

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

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

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:185
void Factorize(I m, std::vector< I > *factors)
Definition: kaldi-math.h:325

## ◆ ComplexFt()

 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.

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) {
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:185
#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.

## ◆ ComplexImExp()

 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.

46  {
47  *a_re = std::cos(x);
48  *a_im = std::sin(x);
49 }

## ◆ ComplexMul()

 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.

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 }

## ◆ ComputeDctMatrix()

 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.

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:98
struct rnnlm::@11::@12 n
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185

## ◆ ComputePca()

 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 616 of file matrix-functions.cc.

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

620  {
621  // Note that some of these matrices may be transposed w.r.t. the
622  // way it's most natural to describe them in math... it's the rows
623  // of X and U that correspond to the (data-points, basis elements).
624  MatrixIndexT N = X.NumRows(), D = X.NumCols();
625  // N = #points, D = feature dim.
626  KALDI_ASSERT(U != NULL && U->NumCols() == D);
627  MatrixIndexT G = U->NumRows(); // # of retained basis elements.
628  KALDI_ASSERT(A == NULL || (A->NumRows() == N && A->NumCols() == G));
629  KALDI_ASSERT(G <= N && G <= D);
630  if (D < N) { // Do conventional PCA.
631  SpMatrix<Real> Msp(D); // Matrix of outer products.
632  Msp.AddMat2(1.0, X, kTrans, 0.0); // M <-- X^T X
633  Matrix<Real> Utmp;
634  Vector<Real> l;
635  if (exact) {
636  Utmp.Resize(D, D);
637  l.Resize(D);
638  //Matrix<Real> M(Msp);
639  //M.DestructiveSvd(&l, &Utmp, NULL);
640  Msp.Eig(&l, &Utmp);
641  } else {
642  Utmp.Resize(D, G);
643  l.Resize(G);
644  Msp.TopEigs(&l, &Utmp);
645  }
646  SortSvd(&l, &Utmp);
647
648  for (MatrixIndexT g = 0; g < G; g++)
649  U->Row(g).CopyColFromMat(Utmp, g);
650  if (print_eigs)
651  KALDI_LOG << (exact ? "" : "Retained ")
652  << "PCA eigenvalues are " << l;
653  if (A != NULL)
654  A->AddMatMat(1.0, X, kNoTrans, *U, kTrans, 0.0);
655  } else { // Do inner-product PCA.
656  SpMatrix<Real> Nsp(N); // Matrix of inner products.
657  Nsp.AddMat2(1.0, X, kNoTrans, 0.0); // M <-- X X^T
658
659  Matrix<Real> Vtmp;
660  Vector<Real> l;
661  if (exact) {
662  Vtmp.Resize(N, N);
663  l.Resize(N);
664  Matrix<Real> Nmat(Nsp);
665  Nmat.DestructiveSvd(&l, &Vtmp, NULL);
666  } else {
667  Vtmp.Resize(N, G);
668  l.Resize(G);
669  Nsp.TopEigs(&l, &Vtmp);
670  }
671
672  MatrixIndexT num_zeroed = 0;
673  for (MatrixIndexT g = 0; g < G; g++) {
674  if (l(g) < 0.0) {
675  KALDI_WARN << "In PCA, setting element " << l(g) << " to zero.";
676  l(g) = 0.0;
677  num_zeroed++;
678  }
679  }
680  SortSvd(&l, &Vtmp); // Make sure zero elements are last, this
681  // is necessary for Orthogonalize() to work properly later.
682
683  Vtmp.Transpose(); // So eigenvalues are the rows.
684
685  for (MatrixIndexT g = 0; g < G; g++) {
686  Real sqrtlg = sqrt(l(g));
687  if (l(g) != 0.0) {
688  U->Row(g).AddMatVec(1.0 / sqrtlg, X, kTrans, Vtmp.Row(g), 0.0);
689  } else {
690  U->Row(g).SetZero();
691  (*U)(g, g) = 1.0; // arbitrary direction. Will later orthogonalize.
692  }
693  if (A != NULL)
694  for (MatrixIndexT n = 0; n < N; n++)
695  (*A)(n, g) = sqrtlg * Vtmp(g, n);
696  }
697  // Now orthogonalize. This is mainly useful in
698  // case there were zero eigenvalues, but we do it
699  // for all of them.
700  U->OrthogonalizeRows();
701  if (print_eigs)
702  KALDI_LOG << "(inner-product) PCA eigenvalues are " << l;
703  }
704 }
int32 MatrixIndexT
Definition: matrix-common.h:98
struct rnnlm::@11::@12 n
#define KALDI_WARN
Definition: kaldi-error.h:150
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
#define KALDI_LOG
Definition: kaldi-error.h:153
void SortSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt, bool sort_on_absolute_value)
Function to ensure that SVD is sorted.

## ◆ CreateEigenvalueMatrix()

 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 2623 of file kaldi-matrix.cc.

Referenced by kaldi::TraceMat(), and kaldi::UnitTestEig().

2624  {
2625  MatrixIndexT n = re.Dim();
2626  KALDI_ASSERT(im.Dim() == n && D->NumRows() == n && D->NumCols() == n);
2627
2628  MatrixIndexT j = 0;
2629  D->SetZero();
2630  while (j < n) {
2631  if (im(j) == 0) { // Real eigenvalue
2632  (*D)(j, j) = re(j);
2633  j++;
2634  } else { // First of a complex pair
2635  KALDI_ASSERT(j+1 < n && ApproxEqual(im(j+1), -im(j))
2636  && ApproxEqual(re(j+1), re(j)));
2638  Real lambda = re(j), mu = im(j);
2639  // create 2x2 block [lambda, mu; -mu, lambda]
2640  (*D)(j, j) = lambda;
2641  (*D)(j, j+1) = mu;
2642  (*D)(j+1, j) = -mu;
2643  (*D)(j+1, j+1) = lambda;
2644  j += 2;
2645  }
2646  }
2647 }
int32 MatrixIndexT
Definition: matrix-common.h:98
struct rnnlm::@11::@12 n
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
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:265

## ◆ RealFft()

 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.

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:98
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:185
#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.

## ◆ RealFftInefficient()

 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.

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:98
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185

 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 :

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 827 of file sp-matrix.cc.

833  {
834  KALDI_ASSERT(Q1.NumRows() == M->NumCols() && P1.NumRows() == M->NumRows() &&
835  G.NumRows() == M->NumRows() && G.NumCols() == M->NumCols() &&
836  M->NumCols() != 0 && Q2.NumRows() == M->NumCols() &&
837  P2.NumRows() == M->NumRows());
838  MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
839  // The following check should not fail as we stipulate P1, P2 and one of Q1
840  // or Q2 must be +ve def and other Q1 or Q2 must be +ve semidef.
841  TpMatrix<Real> LInv(rows);
842  LInv.Cholesky(P1);
843  LInv.Invert(); // Will throw exception if fails.
844  SpMatrix<Real> S(rows);
845  Matrix<Real> LInvFull(LInv);
846  S.AddMat2Sp(1.0, LInvFull, kNoTrans, P2, 0.0); // S := L^{-1} P_2 L^{-T}
847  Matrix<Real> U(rows, rows);
848  Vector<Real> d(rows);
849  S.SymPosSemiDefEig(&d, &U);
850  Matrix<Real> T(rows, rows);
851  T.AddMatMat(1.0, U, kTrans, LInvFull, kNoTrans, 0.0); // T := U^T * L^{-1}
852
853 #ifdef KALDI_PARANOID // checking mainly for errors in the code or math.
854  {
855  SpMatrix<Real> P1Trans(rows);
856  P1Trans.AddMat2Sp(1.0, T, kNoTrans, P1, 0.0);
857  KALDI_ASSERT(P1Trans.IsUnit(0.01));
858  }
859  {
860  SpMatrix<Real> P2Trans(rows);
861  P2Trans.AddMat2Sp(1.0, T, kNoTrans, P2, 0.0);
862  KALDI_ASSERT(P2Trans.IsDiagonal(0.01));
863  }
864 #endif
865
866  Matrix<Real> TInv(T);
867  TInv.Invert();
868  Matrix<Real> Gdash(rows, cols);
869  Gdash.AddMatMat(1.0, T, kNoTrans, G, kNoTrans, 0.0); // G' = T G
870  Matrix<Real> MdashOld(rows, cols);
871  MdashOld.AddMatMat(1.0, TInv, kTrans, *M, kNoTrans, 0.0); // M' = T^{-T} M
872  Matrix<Real> MdashNew(MdashOld);
873  Real objf_impr = 0.0;
874  for (MatrixIndexT n = 0; n < rows; n++) {
875  SpMatrix<Real> Qsum(Q1);
877  SubVector<Real> mdash_n = MdashNew.Row(n);
878  SubVector<Real> gdash_n = Gdash.Row(n);
879
880  Matrix<Real> QsumInv(Qsum);
881  try {
882  QsumInv.Invert();
883  Real old_objf = VecVec(mdash_n, gdash_n)
884  - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
885  mdash_n.AddMatVec(1.0, QsumInv, kNoTrans, gdash_n, 0.0); // m'_n := g'_n * (Q_1 + d_n Q_2)^{-1}
886  Real new_objf = VecVec(mdash_n, gdash_n)
887  - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
888  if (new_objf < old_objf) {
889  if (new_objf < old_objf - 1.0e-05) {
890  KALDI_WARN << "In double quadratic matrix problem: objective "
891  "function decreasing during optimization of " << opts.name
892  << ", " << old_objf << "->" << new_objf << ", change is "
893  << (new_objf - old_objf);
894  KALDI_ERR << "Auxiliary function decreasing."; // Will be caught.
895  } else { // Reset to old value, didn't improve (very close to optimum).
896  MdashNew.Row(n).CopyFromVec(MdashOld.Row(n));
897  }
898  }
899  objf_impr += new_objf - old_objf;
900  }
901  catch (...) {
902  KALDI_WARN << "Matrix inversion or optimization failed during double "
903  "quadratic problem, solving for" << opts.name
904  << ": trying more stable approach.";
905  objf_impr += SolveQuadraticProblem(Qsum, gdash_n, opts, &mdash_n);
906  }
907  }
908  M->AddMatMat(1.0, T, kTrans, MdashNew, kNoTrans, 0.0); // M := T^T M'.
909  return objf_impr;
910 }
double SolveQuadraticProblem(const SpMatrix< double > &H, const VectorBase< double > &g, const SolverOptions &opts, VectorBase< double > *x)
Definition: sp-matrix.cc:635
template double VecSpVec(const VectorBase< double > &v1, const SpMatrix< double > &M, const VectorBase< double > &v2)
int32 MatrixIndexT
Definition: matrix-common.h:98
struct rnnlm::@11::@12 n
#define KALDI_ERR
Definition: kaldi-error.h:147
#define KALDI_WARN
Definition: kaldi-error.h:150
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
Real VecVec(const VectorBase< Real > &a, const VectorBase< Real > &b)
Returns dot product between v1 and v2.
Definition: kaldi-vector.cc:37

 Real SolveQuadraticMatrixProblem ( const SpMatrix< Real > & Q, const MatrixBase< Real > & Y, const SpMatrix< Real > & P, const SolverOptions & opts, MatrixBase< Real > * M )

Maximizes the auxiliary function :

Like a numerically stable version of .

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 729 of file sp-matrix.cc.

733  {
734  KALDI_ASSERT(Q.NumRows() == M->NumCols() &&
735  SigmaInv.NumRows() == M->NumRows() && Y.NumRows() == M->NumRows()
736  && Y.NumCols() == M->NumCols() && M->NumCols() != 0);
737  opts.Check();
738  MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
739  if (Q.IsZero(0.0)) {
740  KALDI_WARN << "Zero quadratic term in quadratic matrix problem for "
741  << opts.name << ": leaving it unchanged.";
742  return 0.0;
743  }
744
745  if (opts.diagonal_precondition) {
746  // We can re-cast the problem with a diagonal preconditioner in the space
747  // of Q (columns of M). Helps to improve the condition of Q.
748  Vector<Real> Q_diag(cols);
749  Q_diag.CopyDiagFromSp(Q);
750  Q_diag.ApplyFloor(std::numeric_limits<Real>::min() * 1.0E+3);
751  Vector<Real> Q_diag_sqrt(Q_diag);
752  Q_diag_sqrt.ApplyPow(0.5);
753  Vector<Real> Q_diag_inv_sqrt(Q_diag_sqrt);
754  Q_diag_inv_sqrt.InvertElements();
755  Matrix<Real> M_scaled(*M);
756  M_scaled.MulColsVec(Q_diag_sqrt);
757  Matrix<Real> Y_scaled(Y);
758  Y_scaled.MulColsVec(Q_diag_inv_sqrt);
759  SpMatrix<Real> Q_scaled(cols);
761  Real ans;
762  SolverOptions new_opts(opts);
763  new_opts.diagonal_precondition = false;
764  ans = SolveQuadraticMatrixProblem(Q_scaled, Y_scaled, SigmaInv,
765  new_opts, &M_scaled);
766  M->CopyFromMat(M_scaled);
767  M->MulColsVec(Q_diag_inv_sqrt);
768  return ans;
769  }
770
771  Matrix<Real> Ybar(Y);
772  if (opts.optimize_delta) {
773  Matrix<Real> Qfull(Q);
774  Ybar.AddMatMat(-1.0, *M, kNoTrans, Qfull, kNoTrans, 1.0);
775  } // Ybar = Y - M Q.
776  Matrix<Real> U(cols, cols);
777  Vector<Real> l(cols);
778  Q.SymPosSemiDefEig(&l, &U); // does svd Q = U L V^T and checks that Q == U L U^T to within a tolerance.
779  // floor l.
780  Real f = std::max<Real>(static_cast<Real>(opts.eps), l.Max() / opts.K);
781  MatrixIndexT nfloored = 0;
782  for (MatrixIndexT i = 0; i < cols; i++) { // floor l.
783  if (l(i) < f) { nfloored++; l(i) = f; }
784  }
785  if (nfloored != 0 && opts.print_debug_output)
786  KALDI_LOG << "Solving matrix problem for " << opts.name
787  << ": floored " << nfloored << " eigenvalues. ";
788  Matrix<Real> tmpDelta(rows, cols);
789  tmpDelta.AddMatMat(1.0, Ybar, kNoTrans, U, kNoTrans, 0.0); // tmpDelta = Ybar * U.
790  l.InvertElements(); KALDI_ASSERT(1.0/l.Max() != 0); // check not infinite. eps should take care of this.
791  tmpDelta.MulColsVec(l); // tmpDelta = Ybar * U * \tilde{L}^{-1}
792
793  Matrix<Real> Delta(rows, cols);
794  Delta.AddMatMat(1.0, tmpDelta, kNoTrans, U, kTrans, 0.0); // Delta = Ybar * U * \tilde{L}^{-1} * U^T
795
796  Real auxf_before, auxf_after;
797  SpMatrix<Real> MQM(rows);
798  Matrix<Real> &SigmaInvY(tmpDelta);
799  { Matrix<Real> SigmaInvFull(SigmaInv); SigmaInvY.AddMatMat(1.0, SigmaInvFull, kNoTrans, Y, kNoTrans, 0.0); }
800  { // get auxf_before. Q(x) = tr(M^T SigmaInv Y) - 0.5 tr(SigmaInv M Q M^T).
801  MQM.AddMat2Sp(1.0, *M, kNoTrans, Q, 0.0);
802  auxf_before = TraceMatMat(*M, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
803  }
804
805  Matrix<Real> Mhat(Delta);
806  if (opts.optimize_delta) Mhat.AddMat(1.0, *M); // Mhat = Delta + M.
807
808  { // get auxf_after.
809  MQM.AddMat2Sp(1.0, Mhat, kNoTrans, Q, 0.0);
810  auxf_after = TraceMatMat(Mhat, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
811  }
812
813  if (auxf_after < auxf_before) {
814  if (auxf_after < auxf_before - 1.0e-10)
815  KALDI_WARN << "Optimizing matrix auxiliary function for "
816  << opts.name << ", auxf decreased "
817  << auxf_before << " to " << auxf_after << ", change is "
818  << (auxf_after-auxf_before);
819  return 0.0;
820  } else {
821  M->CopyFromMat(Mhat);
822  return auxf_after - auxf_before;
823  }
824 }
int32 MatrixIndexT
Definition: matrix-common.h:98
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:150
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:185
Real TraceSpSp(const SpMatrix< Real > &A, const SpMatrix< OtherReal > &B)
Returns tr(A B).
Definition: sp-matrix.cc:367
#define KALDI_LOG
Definition: kaldi-error.h:153

 Real kaldi::SolveQuadraticProblem ( const SpMatrix< Real > & H, const VectorBase< Real > & g, const SolverOptions & opts, VectorBase< Real > * x )

Maximizes the auxiliary function

using a numerically stable method.

Like a numerically stable version of Assumes H positive semidefinite. Returns the objective-function change.

## ◆ SortSvd()

 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 2580 of file kaldi-matrix.cc.

2581  {
2583  MatrixIndexT num_singval = s->Dim();
2584  KALDI_ASSERT(U == NULL || U->NumCols() == num_singval);
2585  KALDI_ASSERT(Vt == NULL || Vt->NumRows() == num_singval);
2586
2587  std::vector<std::pair<Real, MatrixIndexT> > vec(num_singval);
2588  // negative because we want revese order.
2589  for (MatrixIndexT d = 0; d < num_singval; d++) {
2590  Real val = (*s)(d),
2591  sort_val = -(sort_on_absolute_value ? std::abs(val) : val);
2592  vec[d] = std::pair<Real, MatrixIndexT>(sort_val, d);
2593  }
2594  std::sort(vec.begin(), vec.end());
2595  Vector<Real> s_copy(*s);
2596  for (MatrixIndexT d = 0; d < num_singval; d++)
2597  (*s)(d) = s_copy(vec[d].second);
2598  if (U != NULL) {
2599  Matrix<Real> Utmp(*U);
2600  MatrixIndexT dim = Utmp.NumRows();
2601  for (MatrixIndexT d = 0; d < num_singval; d++) {
2602  MatrixIndexT oldidx = vec[d].second;
2603  for (MatrixIndexT e = 0; e < dim; e++)
2604  (*U)(e, d) = Utmp(e, oldidx);
2605  }
2606  }
2607  if (Vt != NULL) {
2608  Matrix<Real> Vttmp(*Vt);
2609  for (MatrixIndexT d = 0; d < num_singval; d++)
2610  (*Vt).Row(d).CopyFromVec(Vttmp.Row(vec[d].second));
2611  }
2612 }
int32 MatrixIndexT
Definition: matrix-common.h:98
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185