qr.cc
Go to the documentation of this file.
1 // matrix/qr.cc
2 
3 // Copyright 2012 Johns Hopkins University (author: Daniel Povey)
4 
5 // See ../../COPYING for clarification regarding multiple authors
6 //
7 // Licensed under the Apache License, Version 2.0 (the "License");
8 // you may not use this file except in compliance with the License.
9 // You may obtain a copy of the License at
10 
11 // http://www.apache.org/licenses/LICENSE-2.0
12 
13 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
14 // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
15 // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
16 // MERCHANTABLITY OR NON-INFRINGEMENT.
17 // See the Apache 2 License for the specific language governing permissions and
18 // limitations under the License.
19 
20 #include <limits>
21 
22 #include "matrix/sp-matrix.h"
23 #include "matrix/kaldi-vector.h"
24 #include "matrix/kaldi-matrix.h"
26 #include "matrix/cblas-wrappers.h"
27 
28 // This file contains an implementation of the Symmetric QR Algorithm
29 // for the symmetric eigenvalue problem. See Golub and Van Loan,
30 // 3rd ed., Algorithm 8.3.3.
31 
32 namespace kaldi {
33 
34 
35 /* This is from Golub and Van Loan 3rd ed., sec. 5.1.3,
36  p210.
37  x is the input of dimenson 'dim', v is the output of dimension
38  dim, and beta is a scalar. Note: we use zero-based
39  not one-based indexing. */
40 /*
41 // We are commenting out the function below ("House") because it's not
42 // needed, but we keep it just to show how we came up with HouseBackward.
43 template<typename Real>
44 void House(MatrixIndexT dim, const Real *x, Real *v, Real *beta) {
45  KALDI_ASSERT(dim > 0);
46  // To avoid overflow, we first compute the max of x_ (or
47  // one if that's zero, and we'll replace "x" by x/max(x_i)
48  // below. The householder vector is anyway invariant to
49  // the magnitude of x. We could actually avoid this extra loop
50  // over x if we wanted to be a bit smarter, but anyway this
51  // doesn't dominate the O(N) performance of the algorithm.
52  Real s; // s is a scale on x.
53  {
54  Real max_x = std::numeric_limits<Real>::min();
55  for (MatrixIndexT i = 0; i < dim; i++)
56  max_x = std::max(max_x, (x[i] < 0 ? -x[i] : x[i]));
57  if (max_x == 0.0) max_x = 1.0;
58  s = 1.0 / max_x;
59  }
60 
61  Real sigma = 0.0;
62  v[0] = 1.0;
63  for (MatrixIndexT i = 1; i < dim; i++) {
64  sigma += (x[i]*s) * (x[i]*s);
65  v[i] = x[i]*s;
66  }
67  if (sigma == 0.0) *beta = 0.0;
68  else {
69  // When we say x1 = x[0], we reference the one-based indexing
70  // in Golub and Van Loan.
71  Real x1 = x[0] * s, mu = std::sqrt(x1*x1 + sigma);
72  if (x1 <= 0) {
73  v[0] = x1 - mu;
74  } else {
75  v[0] = -sigma / (x1 + mu);
76  KALDI_ASSERT(KALDI_ISFINITE(v[dim-1]));
77  }
78  Real v1 = v[0];
79  Real v1sq = v1 * v1;
80  *beta = 2 * v1sq / (sigma + v1sq);
81  Real inv_v1 = 1.0 / v1;
82  if (KALDI_ISINF(inv_v1)) {
83  // can happen if v1 is denormal.
84  KALDI_ASSERT(v1 == v1 && v1 != 0.0);
85  for (MatrixIndexT i = 0; i < dim; i++) v[i] /= v1;
86  } else {
87  cblas_Xscal(dim, inv_v1, v, 1);
88  }
89  if (KALDI_ISNAN(inv_v1)) {
90  KALDI_ERR << "NaN encountered in HouseBackward";
91  }
92  }
93 }
94 */
95 
96 // This is a backward version of the "House" routine above:
97 // backward because it's the last index, not the first index of
98 // the vector that is "special". This is convenient in
99 // the Tridiagonalize routine that uses reversed indexes for
100 // compatibility with the packed lower triangular format.
101 template<typename Real>
102 void HouseBackward(MatrixIndexT dim, const Real *x, Real *v, Real *beta) {
103  KALDI_ASSERT(dim > 0);
104  // To avoid overflow, we first compute the max of x_ (or
105  // one if that's zero, and we'll replace "x" by x/max(x_i)
106  // below. The householder vector is anyway invariant to
107  // the magnitude of x. We could actually avoid this extra loop
108  // over x if we wanted to be a bit smarter, but anyway this
109  // doesn't dominate the O(N) performance of the algorithm.
110  Real s; // s is a scale on x.
111  {
112  Real max_x = std::numeric_limits<Real>::min();
113  for (MatrixIndexT i = 0; i < dim; i++)
114  max_x = std::max(max_x, (x[i] < 0 ? -x[i] : x[i]));
115  s = 1.0 / max_x;
116  }
117  Real sigma = 0.0;
118  v[dim-1] = 1.0;
119  for (MatrixIndexT i = 0; i + 1 < dim; i++) {
120  sigma += (x[i] * s) * (x[i] * s);
121  v[i] = x[i] * s;
122  }
123  KALDI_ASSERT(KALDI_ISFINITE(sigma) &&
124  "Tridiagonalizing matrix that is too large or has NaNs.");
125  if (sigma == 0.0) *beta = 0.0;
126  else {
127  Real x1 = x[dim-1] * s, mu = std::sqrt(x1 * x1 + sigma);
128  if (x1 <= 0) {
129  v[dim-1] = x1 - mu;
130  } else {
131  v[dim-1] = -sigma / (x1 + mu);
132  KALDI_ASSERT(KALDI_ISFINITE(v[dim-1]));
133  }
134  Real v1 = v[dim-1];
135  Real v1sq = v1 * v1;
136  *beta = 2 * v1sq / (sigma + v1sq);
137  Real inv_v1 = 1.0 / v1;
138  if (KALDI_ISINF(inv_v1)) {
139  // can happen if v1 is denormal.
140  KALDI_ASSERT(v1 == v1 && v1 != 0.0);
141  for (MatrixIndexT i = 0; i < dim; i++) v[i] /= v1;
142  } else {
143  cblas_Xscal(dim, inv_v1, v, 1);
144  }
145  if (KALDI_ISNAN(inv_v1)) {
146  KALDI_ERR << "NaN encountered in HouseBackward";
147  }
148  }
149 }
150 
151 
164 template<typename Real>
166  MatrixIndexT n = this->NumRows();
167  KALDI_ASSERT(Q == NULL || (Q->NumRows() == n &&
168  Q->NumCols() == n));
169  if (Q != NULL) Q->SetUnit();
170  Real *data = this->Data();
171  Real *qdata = (Q == NULL ? NULL : Q->Data());
172  MatrixIndexT qstride = (Q == NULL ? 0 : Q->Stride());
173  Vector<Real> tmp_v(n-1), tmp_p(n);
174  Real beta, *v = tmp_v.Data(), *p = tmp_p.Data(), *w = p, *x = p;
175  for (MatrixIndexT k = n-1; k >= 2; k--) {
176  MatrixIndexT ksize = ((k+1)*k)/2;
177  // ksize is the packed size of the lower-triangular matrix of size k,
178  // which is the size of "all rows previous to this one."
179  Real *Arow = data + ksize; // In Golub+Van Loan it was A(k+1:n, k), we
180  // have Arow = A(k, 0:k-1).
181  HouseBackward(k, Arow, v, &beta); // sets v and beta.
182  cblas_Xspmv(k, beta, data, v, 1, 0.0, p, 1); // p = beta * A(0:k-1,0:k-1) v
183  Real minus_half_beta_pv = -0.5 * beta * cblas_Xdot(k, p, 1, v, 1);
184  cblas_Xaxpy(k, minus_half_beta_pv, v, 1, w, 1); // w = p - (beta p^T v/2) v;
185  // this relies on the fact that w and p are the same pointer.
186  // We're doing A(k, k-1) = ||Arow||. It happens that this element
187  // is indexed at ksize + k - 1 in the packed lower-triangular format.
188  data[ksize + k - 1] = std::sqrt(cblas_Xdot(k, Arow, 1, Arow, 1));
189  for (MatrixIndexT i = 0; i + 1 < k; i++)
190  data[ksize + i] = 0; // This is not in Golub and Van Loan but is
191  // necessary if we're not using parts of A to store the Householder
192  // vectors.
193  // We're doing A(0:k-1,0:k-1) -= (v w' + w v')
194  cblas_Xspr2(k, -1.0, v, 1, w, 1, data);
195  if (Q != NULL) { // C.f. Golub, Q is H_1 .. H_n-2... in this
196  // case we apply them in the opposite order so it's H_n-1 .. H_1,
197  // but also Q is transposed so we really have Q = H_1 .. H_n-1.
198  // It's a double negative.
199  // Anyway, we left-multiply Q by each one. The H_n would each be
200  // diag(I + beta v v', I) but we don't ever touch the last dims.
201  // We do (in Matlab notation):
202  // Q(0:k-1,:) = (I - beta v v') * Q, i.e.:
203  // Q(:,0:i-1) += -beta v (v' Q(:,0:k-1)v .. let x = -beta Q(0:k-1,:)^T v.
204  cblas_Xgemv(kTrans, k, n, -beta, qdata, qstride, v, 1, 0.0, x, 1);
205  // now x = -beta Q(:,0:k-1) v.
206  // The next line does: Q(:,0:k-1) += v x'.
207  cblas_Xger(k, n, 1.0, v, 1, x, 1, qdata, qstride);
208  }
209  }
210 }
211 
212 // Instantiate these functions, as it wasn't implemented in sp-matrix.cc
213 // where we instantiated the whole class.
214 template
216 template
218 
220 template<typename Real>
221 inline void Givens(Real a, Real b, Real *c, Real *s) {
222  if (b == 0) {
223  *c = 1;
224  *s = 0;
225  } else {
226  if (std::abs(b) > std::abs(a)) {
227  Real tau = -a / b;
228  *s = 1 / std::sqrt(1 + tau*tau);
229  *c = *s * tau;
230  } else {
231  Real tau = -b / a;
232  *c = 1 / std::sqrt(1 + tau*tau);
233  *s = *c * tau;
234  }
235  }
236 }
237 
238 
239 // Some internal code for the QR algorithm: one "QR step".
240 // This is Golub and Van Loan 3rd ed., Algorithm 8.3.2 "Implicit Symmetric QR step
241 // with Wilkinson shift." A couple of differences: this code is
242 // in zero based arithmetic, and we represent Q transposed from
243 // their Q for memory locality with row-major-indexed matrices.
244 template <typename Real>
246  Real *diag,
247  Real *off_diag,
248  MatrixBase<Real> *Q) {
249  KALDI_ASSERT(n >= 2);
250  // below, "scale" could be any number; we introduce it to keep the
251  // floating point quantities within a good range.
252  Real d = (diag[n-2] - diag[n-1]) / 2.0,
253  t = off_diag[n-2],
254  inv_scale = std::max(std::max(std::abs(d), std::abs(t)),
255  std::numeric_limits<Real>::min()),
256  scale = 1.0 / inv_scale,
257  d_scaled = d * scale,
258  off_diag_n2_scaled = off_diag[n-2] * scale,
259  t2_n_n1_scaled = off_diag_n2_scaled * off_diag_n2_scaled,
260  sgn_d = (d > 0.0 ? 1.0 : -1.0),
261  mu = diag[n-1] - inv_scale * t2_n_n1_scaled /
262  (d_scaled + sgn_d * std::sqrt(d_scaled * d_scaled + t2_n_n1_scaled)),
263  x = diag[0] - mu,
264  z = off_diag[0];
266  Real *Qdata = (Q == NULL ? NULL : Q->Data());
267  MatrixIndexT Qstride = (Q == NULL ? 0 : Q->Stride()),
268  Qcols = (Q == NULL ? 0 : Q->NumCols());
269  for (MatrixIndexT k = 0; k < n-1; k++) {
270  Real c, s;
271  Givens(x, z, &c, &s);
272  // Rotate dimensions k and k+1 with the Givens matrix G, as
273  // T <== G^T T G.
274  // In 2d, a Givens matrix is [ c s; -s c ]. Forget about
275  // the dimension-indexing issues and assume we have a 2x2
276  // symmetric matrix [ p q ; q r ]
277  // We ask our friends at Wolfram Alpha about
278  // { { c, -s}, {s, c} } * { {p, q}, {q, r} } * { { c, s}, {-s, c} }
279  // Interpreting the result as [ p', q' ; q', r ]
280  // p' = c (c p - s q) - s (c q - s r)
281  // q' = s (c p - s q) + c (c q - s r)
282  // r' = s (s p + c q) + c (s q + c r)
283  Real p = diag[k], q = off_diag[k], r = diag[k+1];
284  // p is element k,k; r is element k+1,k+1; q is element k,k+1 or k+1,k.
285  // We'll let the compiler optimize this.
286  diag[k] = c * (c*p - s*q) - s * (c*q - s*r);
287  off_diag[k] = s * (c*p - s*q) + c * (c*q - s*r);
288  diag[k+1] = s * (s*p + c*q) + c * (s*q + c*r);
289 
290  // We also have some other elements to think of that
291  // got rotated in a simpler way: if k>0,
292  // then element (k, k-1) and (k+1, k-1) get rotated. Here,
293  // element k+1, k-1 will be present as z; it's the out-of-band
294  // element that we remembered from last time. This is
295  // on the left as it's the row indexes that differ, so think of
296  // this as being premultiplied by G^T. In fact we're multiplying
297  // T by in some sense the opposite/transpose of the Givens rotation.
298  if (k > 0) { // Note, in rotations, going backward, (x,y) -> ((cx - sy), (sx + cy))
299  Real &elem_k_km1 = off_diag[k-1],
300  elem_kp1_km1 = z; // , tmp = elem_k_km1;
301  elem_k_km1 = c*elem_k_km1 - s*elem_kp1_km1;
302  // The next line will set elem_kp1_km1 to zero and we'll never access this
303  // value, so we comment it out.
304  // elem_kp1_km1 = s*tmp + c*elem_kp1_km1;
305  }
306  if (Q != NULL)
307  cblas_Xrot(Qcols, Qdata + k*Qstride, 1,
308  Qdata + (k+1)*Qstride, 1, c, -s);
309  if (k < n-2) {
310  // Next is the elements (k+2, k) and (k+2, k-1), to be rotated, again
311  // backwards.
312  Real &elem_kp2_k = z,
313  &elem_kp2_kp1 = off_diag[k+1];
314  // Note: elem_kp2_k == z would start off as zero because it's
315  // two off the diagonal, and not been touched yet. Therefore
316  // we eliminate it in expressions below, commenting it out.
317  // If we didn't do this we should set it to zero first.
318  elem_kp2_k = - s * elem_kp2_kp1; // + c*elem_kp2_k
319  elem_kp2_kp1 = c * elem_kp2_kp1; // + s*elem_kp2_k (original value).
320  // The next part is from the algorithm they describe: x = t_{k+1,k}
321  x = off_diag[k];
322  }
323  }
324 }
325 
326 
327 // Internal code for the QR algorithm, where the diagonal
328 // and off-diagonal of the symmetric matrix are represented as
329 // vectors of length n and n-1.
330 template <typename Real>
332  Real *diag,
333  Real *off_diag,
334  MatrixBase<Real> *Q) {
335  KALDI_ASSERT(Q == NULL || Q->NumCols() == n); // We may
336  // later relax the condition that Q->NumCols() == n.
337 
338  MatrixIndexT counter = 0, max_iters = 500 + 4*n, // Should never take this many iters.
339  large_iters = 100 + 2*n;
340  Real epsilon = (pow(2.0, sizeof(Real) == 4 ? -23.0 : -52.0));
341 
342  for (; counter < max_iters; counter++) { // this takes the place of "until
343  // q=n"... we'll break out of the
344  // loop when we converge.
345  if (counter == large_iters ||
346  (counter > large_iters && (counter - large_iters) % 50 == 0)) {
347  KALDI_WARN << "Took " << counter
348  << " iterations in QR (dim is " << n << "), doubling epsilon.";
349  SubVector<Real> d(diag, n), o(off_diag, n-1);
350  KALDI_WARN << "Diag, off-diag are " << d << " and " << o;
351  epsilon *= 2.0;
352  }
353  for (MatrixIndexT i = 0; i+1 < n; i++) {
354  if (std::abs(off_diag[i]) <= epsilon *
355  (std::abs(diag[i]) + std::abs(diag[i+1])))
356  off_diag[i] = 0.0;
357  }
358  // The next code works out p, q, and npq which is n - p - q.
359  // For the definitions of q and p, see Golub and Van Loan; we
360  // partition the n dims into pieces of size (p, n-p-q, q) where
361  // the part of size q is diagonal and the part of size n-p-p is
362  // "unreduced", i.e. has no zero off-diagonal elements.
363  MatrixIndexT q = 0;
364  // Note: below, "n-q < 2" should more clearly be "n-2-q < 0", but that
365  // causes problems if MatrixIndexT is unsigned.
366  while (q < n && (n-q < 2 || off_diag[n-2-q] == 0.0))
367  q++;
368  if (q == n) break; // we're done. It's diagonal.
369  KALDI_ASSERT(n - q >= 2);
370  MatrixIndexT npq = 2; // Value of n - p - q, where n - p - q must be
371  // unreduced. This is the size of "middle" band of elements. If q != n,
372  // we must have hit a nonzero off-diag element, so the size of this
373  // band must be at least two.
374  while (npq + q < n && (n-q-npq-1 < 0 || off_diag[n-q-npq-1] != 0.0))
375  npq++;
376  MatrixIndexT p = n - q - npq;
377  { // Checks.
378  for (MatrixIndexT i = 0; i+1 < npq; i++)
379  KALDI_ASSERT(off_diag[p + i] != 0.0);
380  for (MatrixIndexT i = 0; i+1 < q; i++)
381  KALDI_ASSERT(off_diag[p + npq - 1 + i] == 0.0);
382  if (p > 1) // Something must have stopped npq from growing further..
383  KALDI_ASSERT(off_diag[p-1] == 0.0); // so last off-diag elem in
384  // group of size p must be zero.
385  }
386 
387  if (Q != NULL) {
388  // Do one QR step on the middle part of Q only.
389  // Qpart will be a subset of the rows of Q.
390  SubMatrix<Real> Qpart(*Q, p, npq, 0, Q->NumCols());
391  QrStep(npq, diag + p, off_diag + p, &Qpart);
392  } else {
393  QrStep(npq, diag + p, off_diag + p,
394  static_cast<MatrixBase<Real>*>(NULL));
395  }
396  }
397  if (counter == max_iters) {
398  KALDI_WARN << "Failure to converge in QR algorithm. "
399  << "Exiting with partial output.";
400  }
401 }
402 
403 
408 template <typename Real>
410  KALDI_ASSERT(this->IsTridiagonal());
411  // We envisage that Q would be square but we don't check for this,
412  // as there are situations where you might not want this.
413  KALDI_ASSERT(Q == NULL || Q->NumRows() == this->NumRows());
414  // Note: the first couple of lines of the algorithm they give would be done
415  // outside of this function, by calling Tridiagonalize().
416 
417  MatrixIndexT n = this->NumRows();
418  Vector<Real> diag(n), off_diag(n-1);
419  for (MatrixIndexT i = 0; i < n; i++) {
420  diag(i) = (*this)(i, i);
421  if (i > 0) off_diag(i-1) = (*this)(i, i-1);
422  }
423  QrInternal(n, diag.Data(), off_diag.Data(), Q);
424  // Now set *this to the value represented by diag and off_diag.
425  this->SetZero();
426  for (MatrixIndexT i = 0; i < n; i++) {
427  (*this)(i, i) = diag(i);
428  if (i > 0) (*this)(i, i-1) = off_diag(i-1);
429  }
430 }
431 
432 template<typename Real>
434  MatrixIndexT dim = this->NumRows();
435  KALDI_ASSERT(s->Dim() == dim);
436  KALDI_ASSERT(P == NULL || (P->NumRows() == dim && P->NumCols() == dim));
437 
438  SpMatrix<Real> A(*this); // Copy *this, since the tridiagonalization
439  // and QR decomposition are destructive.
440  // Note: for efficiency of memory access, the tridiagonalization
441  // algorithm makes the *rows* of P the eigenvectors, not the columns.
442  // We'll transpose P before we exit.
443  // Also note: P may be null if you don't want the eigenvectors. This
444  // will make this function more efficient.
445 
446  A.Tridiagonalize(P); // Tridiagonalizes.
447  A.Qr(P); // Diagonalizes.
448  if(P) P->Transpose();
449  s->CopyDiagFromPacked(A);
450 }
451 
452 
453 template<typename Real>
455  MatrixIndexT lanczos_dim) const {
456  const SpMatrix<Real> &S(*this); // call this "S" for easy notation.
457  MatrixIndexT eig_dim = s->Dim(); // Space of dim we want to retain.
458  if (lanczos_dim <= 0)
459  lanczos_dim = std::max(eig_dim + 50, eig_dim + eig_dim/2);
460  MatrixIndexT dim = this->NumRows();
461  if (lanczos_dim >= dim) {
462  // There would be no speed advantage in using this method, so just
463  // use the regular approach.
464  Vector<Real> s_tmp(dim);
465  Matrix<Real> P_tmp(dim, dim);
466  this->Eig(&s_tmp, &P_tmp);
467  SortSvd(&s_tmp, &P_tmp);
468  s->CopyFromVec(s_tmp.Range(0, eig_dim));
469  P->CopyFromMat(P_tmp.Range(0, dim, 0, eig_dim));
470  return;
471  }
472  KALDI_ASSERT(eig_dim <= dim && eig_dim > 0);
473  KALDI_ASSERT(P->NumRows() == dim && P->NumCols() == eig_dim); // each column
474  // is one eigenvector.
475 
476  Matrix<Real> Q(lanczos_dim, dim); // The rows of Q will be the
477  // orthogonal vectors of the Krylov subspace.
478 
479  SpMatrix<Real> T(lanczos_dim); // This will be equal to Q S Q^T,
480  // i.e. *this projected into the Krylov subspace. Note: only the
481  // diagonal and off-diagonal fo T are nonzero, i.e. it's tridiagonal,
482  // but we don't have access to the low-level algorithms that work
483  // on that type of matrix (since we want to use ATLAS). So we just
484  // do normal SVD, on a full matrix; it won't typically dominate.
485 
486  Q.Row(0).SetRandn();
487  Q.Row(0).Scale(1.0 / Q.Row(0).Norm(2));
488  for (MatrixIndexT d = 0; d < lanczos_dim; d++) {
489  Vector<Real> r(dim);
490  r.AddSpVec(1.0, S, Q.Row(d), 0.0);
491  // r = S * q_d
492  MatrixIndexT counter = 0;
493  Real end_prod;
494  while (1) { // Normally we'll do this loop only once:
495  // we repeat to handle cases where r gets very much smaller
496  // and we want to orthogonalize again.
497  // We do "full orthogonalization" to preserve stability,
498  // even though this is usually a waste of time.
499  Real start_prod = VecVec(r, r);
500  for (SignedMatrixIndexT e = d; e >= 0; e--) { // e must be signed!
501  SubVector<Real> q_e(Q, e);
502  Real prod = VecVec(r, q_e);
503  if (counter == 0 && static_cast<MatrixIndexT>(e) + 1 >= d) // Keep T tridiagonal, which
504  T(d, e) = prod; // mathematically speaking, it is.
505  r.AddVec(-prod, q_e); // Subtract component in q_e.
506  }
507  if (d+1 == lanczos_dim) break;
508  end_prod = VecVec(r, r);
509  if (end_prod <= 0.1 * start_prod) {
510  // also handles case where both are 0.
511  // We're not confident any more that it's completely
512  // orthogonal to the rest so we want to re-do.
513  if (end_prod == 0.0)
514  r.SetRandn(); // "Restarting".
515  counter++;
516  if (counter > 100)
517  KALDI_ERR << "Loop detected in Lanczos iteration.";
518  } else {
519  break;
520  }
521  }
522  if (d+1 != lanczos_dim) {
523  // OK, at this point we're satisfied that r is orthogonal
524  // to all previous rows.
525  KALDI_ASSERT(end_prod != 0.0); // should have looped.
526  r.Scale(1.0 / std::sqrt(end_prod)); // make it unit.
527  Q.Row(d+1).CopyFromVec(r);
528  }
529  }
530 
531  Matrix<Real> R(lanczos_dim, lanczos_dim);
532  R.SetUnit();
533  T.Qr(&R); // Diagonalizes T.
534  Vector<Real> s_tmp(lanczos_dim);
535  s_tmp.CopyDiagFromSp(T);
536 
537  // Now T = R * diag(s_tmp) * R^T.
538  // The next call sorts the elements of s from greatest to least absolute value,
539  // and moves around the rows of R in the corresponding way. This picks out
540  // the largest (absolute) eigenvalues.
541  SortSvd(&s_tmp, static_cast<Matrix<Real>*>(NULL), &R);
542  // Keep only the initial rows of R, those corresponding to greatest (absolute)
543  // eigenvalues.
544  SubMatrix<Real> Rsub(R, 0, eig_dim, 0, lanczos_dim);
545  SubVector<Real> s_sub(s_tmp, 0, eig_dim);
546  s->CopyFromVec(s_sub);
547 
548  // For working out what to do now, just assume the other eigenvalues were
549  // zero. This is just for purposes of knowing how to get the result, and
550  // not getting things wrongly transposed.
551  // We have T = Rsub^T * diag(s_sub) * Rsub.
552  // Now, T = Q S Q^T, with Q orthogonal, so S = Q^T T Q = Q^T Rsub^T * diag(s) * Rsub * Q.
553  // The output is P and we want S = P * diag(s) * P^T, so we need P = Q^T Rsub^T.
554  P->AddMatMat(1.0, Q, kTrans, Rsub, kTrans, 0.0);
555 }
556 
557 
558 // Instantiate the templates for Eig and TopEig.
559 template
561 template
563 
564 template
566 template
568 
569 // Someone had a problem with the Intel compiler with -O3, with Qr not being
570 // defined for some strange reason (should automatically happen when
571 // we instantiate Eig and TopEigs), so we explicitly instantiate it here.
572 template
574 template
576 
577 
578 
579 }
580 // namespace kaldi
This code computes Goodness of Pronunciation (GOP) and extracts phone-level pronunciation feature for...
Definition: chain.dox:20
void Givens(Real a, Real b, Real *c, Real *s)
Create Givens rotations, as in Golub and Van Loan 3rd ed., page 216.
Definition: qr.cc:221
Packed symetric matrix class.
Definition: matrix-common.h:62
void QrInternal(MatrixIndexT n, Real *diag, Real *off_diag, MatrixBase< Real > *Q)
Definition: qr.cc:331
void HouseBackward(MatrixIndexT dim, const Real *x, Real *v, Real *beta)
Definition: qr.cc:102
void cblas_Xspr2(MatrixIndexT dim, float alpha, const float *Xdata, MatrixIndexT incX, const float *Ydata, MatrixIndexT incY, float *Adata)
void Qr(MatrixBase< Real > *Q)
The symmetric QR algorithm.
Definition: qr.cc:409
MatrixIndexT NumCols() const
Returns number of columns (or zero for empty matrix).
Definition: kaldi-matrix.h:67
Base class which provides matrix operations not involving resizing or allocation. ...
Definition: kaldi-matrix.h:49
const Real * Data() const
Gives pointer to raw data (const).
Definition: kaldi-matrix.h:79
#define KALDI_ISINF
Definition: kaldi-math.h:73
#define KALDI_ISFINITE(x)
Definition: kaldi-math.h:74
void Eig(VectorBase< Real > *s, MatrixBase< Real > *P=NULL) const
Solves the symmetric eigenvalue problem: at end we should have (*this) = P * diag(s) * P^T...
Definition: qr.cc:433
A class for storing matrices.
Definition: kaldi-matrix.h:823
void AddSpVec(const Real alpha, const SpMatrix< Real > &M, const VectorBase< Real > &v, const Real beta)
Add symmetric positive definite matrix times vector: this <– beta*this + alpha*M*v.
void CopyFromMat(const MatrixBase< OtherReal > &M, MatrixTransposeType trans=kNoTrans)
Copy given matrix. (no resize is done).
void SetUnit()
Sets to zero, except ones along diagonal [for non-square matrices too].
void Tridiagonalize(MatrixBase< Real > *Q)
Tridiagonalize the matrix with an orthogonal transformation.
Definition: qr.cc:165
float cblas_Xdot(const int N, const float *const X, const int incX, const float *const Y, const int incY)
void CopyFromVec(const VectorBase< Real > &v)
Copy data from another vector (must match own size).
void CopyDiagFromPacked(const PackedMatrix< Real > &M)
Extracts the diagonal of a packed matrix M; works for Sp or Tp.
MatrixIndexT Stride() const
Stride (distance in memory between each row). Will be >= NumCols.
Definition: kaldi-matrix.h:70
int32 MatrixIndexT
Definition: matrix-common.h:98
void cblas_Xscal(const int N, const float alpha, float *data, const int inc)
struct rnnlm::@11::@12 n
void Transpose()
Transpose the matrix.
void AddMatMat(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const Real beta)
#define KALDI_ERR
Definition: kaldi-error.h:147
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 corre...
Definition: qr.cc:454
#define KALDI_WARN
Definition: kaldi-error.h:150
Real * Data()
Returns a pointer to the start of the vector&#39;s data.
Definition: kaldi-vector.h:70
MatrixIndexT Dim() const
Returns the dimension of the vector.
Definition: kaldi-vector.h:64
void Scale(Real alpha)
Multiplies all elements by this constant.
void SetRandn()
Set vector to random normally-distributed noise.
A class representing a vector.
Definition: kaldi-vector.h:406
void cblas_Xaxpy(const int N, const float alpha, const float *X, const int incX, float *Y, const int incY)
#define KALDI_ISNAN
Definition: kaldi-math.h:72
int32 SignedMatrixIndexT
Definition: matrix-common.h:99
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
MatrixIndexT NumRows() const
Returns number of rows (or zero for empty matrix).
Definition: kaldi-matrix.h:64
void QrStep(MatrixIndexT n, Real *diag, Real *off_diag, MatrixBase< Real > *Q)
Definition: qr.cc:245
void cblas_Xgemv(MatrixTransposeType trans, MatrixIndexT num_rows, MatrixIndexT num_cols, float alpha, const float *Mdata, MatrixIndexT stride, const float *xdata, MatrixIndexT incX, float beta, float *ydata, MatrixIndexT incY)
SubMatrix< Real > Range(const MatrixIndexT row_offset, const MatrixIndexT num_rows, const MatrixIndexT col_offset, const MatrixIndexT num_cols) const
Return a sub-part of matrix.
Definition: kaldi-matrix.h:202
void cblas_Xspmv(const float alpha, const int num_rows, const float *Mdata, const float *v, const int v_inc, const float beta, float *y, const int y_inc)
Provides a vector abstraction class.
Definition: kaldi-vector.h:41
void cblas_Xger(MatrixIndexT num_rows, MatrixIndexT num_cols, float alpha, const float *xdata, MatrixIndexT incX, const float *ydata, MatrixIndexT incY, float *Mdata, MatrixIndexT stride)
Real VecVec(const VectorBase< Real > &a, const VectorBase< Real > &b)
Returns dot product between v1 and v2.
Definition: kaldi-vector.cc:37
void AddVec(const Real alpha, const VectorBase< OtherReal > &v)
Add vector : *this = *this + alpha * rv (with casting between floats and doubles) ...
void CopyDiagFromSp(const SpMatrix< Real > &M)
Extracts the diagonal of a symmetric matrix.
Definition: kaldi-vector.h:309
Sub-matrix representation.
Definition: kaldi-matrix.h:988
Represents a non-allocating general vector which can be defined as a sub-vector of higher-level vecto...
Definition: kaldi-vector.h:501
void cblas_Xrot(const int N, float *X, const int incX, float *Y, const int incY, const float c, const float s)
void SortSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt, bool sort_on_absolute_value)
Function to ensure that SVD is sorted.
SubVector< Real > Range(const MatrixIndexT o, const MatrixIndexT l)
Returns a sub-vector of a vector (a range of elements).
Definition: kaldi-vector.h:94