matrix-lib-test.cc
Go to the documentation of this file.
1 // matrix/matrix-lib-test.cc
2 
3 // Copyright 2009-2012 Microsoft Corporation; Mohit Agarwal; Lukas Burget;
4 // Ondrej Glembek; Saarland University (Author: Arnab Ghoshal);
5 // Go Vivace Inc.; Yanmin Qian; Jan Silovsky;
6 // Johns Hopkins University (Author: Daniel Povey);
7 // Haihua Xu; Wei Shi
8 // 2015 Guoguo Chen
9 // 2017 Daniel Galvez
10 
11 // See ../../COPYING for clarification regarding multiple authors
12 //
13 // Licensed under the Apache License, Version 2.0 (the "License");
14 // you may not use this file except in compliance with the License.
15 // You may obtain a copy of the License at
16 
17 // http://www.apache.org/licenses/LICENSE-2.0
18 
19 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
20 // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
21 // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
22 // MERCHANTABLITY OR NON-INFRINGEMENT.
23 // See the Apache 2 License for the specific language governing permissions and
24 // limitations under the License.
25 
26 #include "matrix/matrix-lib.h"
27 #include "util/stl-utils.h"
28 #include <numeric>
29 #include <time.h> // This is only needed for UnitTestSvdSpeed, you can
30 // comment it (and that function) out if it causes problems.
31 #include <matrix/cblas-wrappers.h>
32 
33 namespace kaldi {
34 
35 template<typename Real>
37  MatrixIndexT dim2 = dim + (Rand() % 3); // slightly higher-dim.
38  // generate random (non-singular) matrix
39  Matrix<Real> tmp(dim, dim2);
40  while (1) {
41  tmp.SetRandn();
42  if (tmp.Cond() < 100) break;
43  KALDI_LOG << "Condition number of random matrix large "
44  << static_cast<float>(tmp.Cond())
45  << ", trying again (this is normal)";
46  }
47  // tmp * tmp^T will give positive definite matrix
48  matrix->AddMat2(1.0, tmp, kNoTrans, 0.0);
49 
50  // Checks that the matrix is indeed pos-def
51  TpMatrix<Real> sqrt(dim);
52  sqrt.Cholesky(*matrix);
53 }
54 
55 
56 template<typename Real> static void InitRandNonsingular(MatrixBase<Real> *M) {
57 start:
58  for (MatrixIndexT i = 0;i < M->NumRows();i++)
59  for (MatrixIndexT j = 0;j < M->NumCols();j++)
60  (*M)(i, j) = RandGauss();
61  if (M->NumRows() != 0 && M->Cond() > 100) {
62  printf("Condition number of random matrix large %f, trying again (this is normal)\n",
63  (float) M->Cond());
64  goto start;
65  }
66 }
67 
68 
69 template<typename Real> static void InitRandNonsingular(SpMatrix<Real> *M) {
70 start:
71  for (MatrixIndexT i = 0;i < M->NumRows();i++)
72  for (MatrixIndexT j = 0;j<=i;j++)
73  (*M)(i, j) = RandGauss();
74  if (M->NumRows() != 0 && M->Cond() > 100)
75  goto start;
76 }
77 
78 /*
79  HERE(2)
80 template<typename Real> static void AssertEqual(const Matrix<Real> &A,
81  const Matrix<Real> &B,
82  float tol = 0.001) {
83  KALDI_ASSERT(A.NumRows() == B.NumRows()&&A.NumCols() == B.NumCols());
84  for (MatrixIndexT i = 0;i < A.NumRows();i++)
85  for (MatrixIndexT j = 0;j < A.NumCols();j++) {
86  KALDI_ASSERT(std::abs(A(i, j)-B(i, j)) < tol*std::max(1.0, (double) (std::abs(A(i, j))+std::abs(B(i, j)))));
87  }
88 }
89 
90 template<typename Real> static void AssertEqual(const SpMatrix<Real> &A,
91  const SpMatrix<Real> &B,
92  float tol = 0.001) {
93  KALDI_ASSERT(A.NumRows() == B.NumRows()&&A.NumCols() == B.NumCols());
94  for (MatrixIndexT i = 0;i < A.NumRows();i++)
95  for (MatrixIndexT j = 0;j<=i;j++)
96  KALDI_ASSERT(std::abs(A(i, j)-B(i, j)) < tol*std::max(1.0, (double) (std::abs(A(i, j))+std::abs(B(i, j)))));
97 }
98 */
99 /*
100  HERE:
101 template<typename Real>
102 static bool ApproxEqual(const SpMatrix<Real> &A,
103  const SpMatrix<Real> &B, Real tol = 0.001) {
104  KALDI_ASSERT(A.NumRows() == B.NumRows());
105  SpMatrix<Real> diff(A);
106  diff.AddSp(1.0, B);
107  Real a = std::max(A.Max(), -A.Min()), b = std::max(B.Max(), -B.Min),
108  d = std::max(diff.Max(), -diff.Min());
109  return (d <= tol * std::max(a, b));
110 }
111 */
112 
113 /* was:
114  template<typename Real>
115  bool ApproxEqual(SpMatrix<Real> &A, SpMatrix<Real> &B, float tol = 0.001) {
116  KALDI_ASSERT(A.NumRows() == B.NumRows()&&A.NumCols() == B.NumCols());
117  for (MatrixIndexT i = 0;i < A.NumRows();i++)
118  for (MatrixIndexT j = 0;j<=i;j++)
119  if (std::abs(A(i, j)-B(i, j)) > tol*std::max(1.0, (double) (std::abs(A(i, j))+std::abs(B(i, j))))) return false;
120  return true;
121  }
122 */
123 
124 /*
125  HERE
126 template<typename Real> static void AssertEqual(Vector<Real> &A, Vector<Real> &B, float tol = 0.001) {
127  KALDI_ASSERT(A.Dim() == B.Dim());
128  for (MatrixIndexT i = 0;i < A.Dim();i++)
129  KALDI_ASSERT(std::abs(A(i)-B(i)) < tol);
130 }
131 
132 template<typename Real> static bool ApproxEqual(Vector<Real> &A, Vector<Real> &B, float tol = 0.001) {
133  KALDI_ASSERT(A.Dim() == B.Dim());
134  for (MatrixIndexT i = 0;i < A.Dim();i++)
135  if (std::abs(A(i)-B(i)) > tol) return false;
136  return true;
137 }
138 */
139 
140 template<typename Real> static void CholeskyUnitTestTr() {
141  for (MatrixIndexT i = 0; i < 5; i++) {
142  MatrixIndexT dimM = 2 + Rand() % 10;
143  Matrix<Real> M(dimM, dimM);
145  SpMatrix<Real> S(dimM);
146  S.AddMat2(1.0, M, kNoTrans, 0.0);
147  TpMatrix<Real> C(dimM);
148  C.Cholesky(S);
149  Matrix<Real> CM(C);
150  TpMatrix<Real> Cinv(C);
151  Cinv.Invert();
152  {
153  Matrix<Real> A(C), B(Cinv), AB(dimM, dimM);
154  AB.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0);
155  KALDI_ASSERT(AB.IsUnit());
156  }
157  SpMatrix<Real> S2(dimM);
158  S2.AddMat2(1.0, CM, kNoTrans, 0.0);
159  AssertEqual(S, S2);
160  C.Invert();
161  Matrix<Real> CM2(C);
162  CM2.Invert();
163  SpMatrix<Real> S3(dimM);
164  S3.AddMat2(1.0, CM2, kNoTrans, 0.0);
165  AssertEqual(S, S3);
166  }
167 }
168 
169 template<typename Real> static void SlowMatMul() {
170  MatrixIndexT N = 1000;
171  Matrix<Real> M(N, N), P(N, N), Q(N, N);
172  for (MatrixIndexT i = 0; i < 10000; i++) {
173  Q.AddMatMat(1.0, M, kNoTrans, P, kNoTrans, 0.0);
174  }
175 }
176 
177 template<typename Real> static void UnitTestAddToDiagMatrix() {
178  for (int p = 0; p < 2; p++) {
179  MatrixIndexT dimM = 10 + Rand() % 2, dimN = 1 + Rand() % 5;
180  Matrix<Real> M(dimM, dimN), Mcopy(M);
181  BaseFloat alpha = 0.35;
182  M.AddToDiag(alpha);
183  for (MatrixIndexT i = 0; i < dimM && i < dimN; i++)
184  Mcopy(i, i) += alpha;
185  AssertEqual(M, Mcopy);
186  }
187 }
188 
189 
190 template<typename Real> static void UnitTestAddDiagVecMat() {
191  for (int p = 0; p < 2; p++) {
192  MatrixIndexT dimM = 100 + Rand() % 255, dimN = 100 + Rand() % 255;
193  Real alpha = 0.43243, beta = 1.423;
194  Matrix<Real> M(dimM, dimN), N(dimM, dimN);
195  M.SetRandn();
196  N.SetRandn();
197  MatrixTransposeType trans = (p % 2 == 0 ? kNoTrans : kTrans);
198  if (trans == kTrans)
199  N.Transpose();
200 
201  Vector<Real> V(dimM);
202  V.SetRandn();
203 
204  Matrix<Real> Mcheck(M);
205  for (int32 r = 0; r < dimM; r++) {
206  SubVector<Real> Mcheckrow(Mcheck, r);
207  Vector<Real> Nrow(dimN);
208  if (trans == kTrans) Nrow.CopyColFromMat(N, r);
209  else Nrow.CopyFromVec(N.Row(r));
210  Mcheckrow.Scale(beta);
211  Mcheckrow.AddVec(alpha * V(r), Nrow);
212  }
213 
214  M.AddDiagVecMat(alpha, V, N, trans, beta);
215  AssertEqual(M, Mcheck);
216  KALDI_ASSERT(M.Sum() != 0.0);
217  }
218 }
219 
220 template<typename Real> static void UnitTestAddMatDiagVec() {
221  // M <- alpha * N[^T] * diag(v) + beta * M
222  for (int p = 0; p < 2; p++) {
223  MatrixIndexT dimM = 100 + Rand() % 255, dimN = 100 + Rand() % 255;
224  Real alpha = 0.43243, beta = 1.423;
225 
226  Matrix<Real> M(dimM, dimN), N(dimM, dimN), buf(dimM, dimN);
227  M.SetRandn();
228  N.SetRandn();
229  buf.CopyFromMat(N);
230  MatrixTransposeType trans = (p % 2 == 0 ? kNoTrans : kTrans);
231  if (trans == kTrans)
232  N.Transpose();
233 
234  Vector<Real> V(dimN);
235  V.SetRandn();
236 
237  Matrix<Real> Mcheck(M);
238  Mcheck.Scale(beta);
239  buf.MulColsVec(V);
240  Mcheck.AddMat(alpha, buf, kNoTrans);
241 
242  M.AddMatDiagVec(alpha, N, trans, V, beta);
243  AssertEqual(M, Mcheck);
244  KALDI_ASSERT(M.Sum() != 0.0);
245  }
246 }
247 
248 template<typename Real> static void UnitTestAddMatMatElements() {
249  // M <- alpha *(A .* B) + beta * M
250  MatrixIndexT dimM = 100 + Rand() % 255, dimN = 100 + Rand() % 255;
251  Real alpha = 0.43243, beta = 1.423;
252  Matrix<Real> M(dimM, dimN), A(dimM, dimN), B(dimM, dimN), buf(dimM, dimN);
253  M.SetRandn();
254  A.SetRandn();
255  B.SetRandn();
256 
257  Matrix<Real> Mcheck(M);
258  buf.CopyFromMat(A); buf.MulElements(B);
259  Mcheck.Scale(beta); Mcheck.AddMat(alpha, buf, kNoTrans);
260 
261  M.AddMatMatElements(alpha, A, B, beta);
262  AssertEqual(M, Mcheck);
263  KALDI_ASSERT(M.Sum() != 0.0);
264 }
265 
266 template<typename Real> static void UnitTestAddSp() {
267  for (MatrixIndexT i = 0;i< 10;i++) {
268  MatrixIndexT dimM = 10+Rand()%10;
269  SpMatrix<Real> S(dimM);
270  S.SetRandn();
271  Matrix<Real> M(S), N(S);
272  N.AddSp(2.0, S);
273  M.Scale(3.0);
274  AssertEqual(M, N);
275  }
276 }
277 
278 template<typename Real, typename OtherReal>
279 static void UnitTestSpAddDiagVec() {
280  for (MatrixIndexT i = 0;i< 10;i++) {
281  BaseFloat alpha = (i<5 ? 1.0 : 0.5);
282  MatrixIndexT dimM = 10+Rand()%10;
283  SpMatrix<Real> S(dimM);
284  S.SetRandn();
285  SpMatrix<Real> T(S);
286  Vector<OtherReal> v(dimM);
287  v.SetRandn();
288  S.AddDiagVec(alpha, v);
289  for (MatrixIndexT i = 0; i < dimM; i++)
290  T(i, i) += alpha * v(i);
291  AssertEqual(S, T);
292  }
293 }
294 
295 
296 template<typename Real>
297 static void UnitTestSpAddVecVec() {
298  for (MatrixIndexT i = 0;i< 10;i++) {
299  BaseFloat alpha = (i<5 ? 1.0 : 0.5);
300  MatrixIndexT dimM = 10+Rand()%10;
301  SpMatrix<Real> S(dimM);
302  S.SetRandn();
303  Matrix<Real> T(S);
304 
305  Vector<Real> v(dimM), w(dimM);
306  v.SetRandn();
307  w.SetRandn();
308  S.AddVecVec(alpha, v, w);
309  T.AddVecVec(alpha, v, w);
310  T.AddVecVec(alpha, w, v);
311  Matrix<Real> U(S);
312  AssertEqual(U, T);
313  }
314 }
315 
316 
317 template<typename Real> static void UnitTestCopyRowsAndCols() {
318  // Test other mode of CopyRowsFromVec, and CopyColsFromVec,
319  // where vector is duplicated.
320  for (MatrixIndexT i = 0; i < 30; i++) {
321  MatrixIndexT dimM = 1 + Rand() % 5, dimN = 1 + Rand() % 5;
322  Vector<float> w(dimN); // test cross-type version of
323  // CopyRowsFromVec.
324  Vector<Real> v(dimM);
325  Matrix<Real> M(dimM, dimN), N(dimM, dimN);
326  v.SetRandn();
327  w.SetRandn();
328  M.CopyColsFromVec(v);
329  N.CopyRowsFromVec(w);
330  for (MatrixIndexT r = 0; r < dimM; r++) {
331  for (MatrixIndexT c = 0; c < dimN; c++) {
332  KALDI_ASSERT(M(r, c) == v(r));
333  KALDI_ASSERT(N(r, c) == w(c));
334  }
335  }
336  }
337 }
338 
339 template<typename Real> static void UnitTestSpliceRows() {
340 
341  for (MatrixIndexT i = 0;i< 10;i++) {
342  MatrixIndexT dimM = 10+Rand()%10, dimN = 10+Rand()%10;
343 
344  Vector<Real> V(dimM*dimN), V10(dimM*dimN);
345  Vector<Real> Vs(std::min(dimM, dimN)), Vs10(std::min(dimM, dimN));
346  V.SetRandn();
347  Matrix<Real> M(dimM, dimN);
348  M.CopyRowsFromVec(V);
349  V10.CopyRowsFromMat(M);
350  AssertEqual(V, V10);
351 
352  for (MatrixIndexT i = 0;i < dimM;i++)
353  for (MatrixIndexT j = 0;j < dimN;j++)
354  KALDI_ASSERT(M(i, j) == V(i*dimN + j));
355 
356  {
357  Vector<Real> V2(dimM), V3(dimM);
358  V2.SetRandn();
359  MatrixIndexT x;
360  M.CopyColFromVec(V2, x = (Rand() % dimN));
361  V3.CopyColFromMat(M, x);
362  AssertEqual(V2, V3);
363  }
364 
365  {
366  Vector<Real> V2(dimN), V3(dimN);
367  V2.SetRandn();
368  MatrixIndexT x;
369  M.CopyRowFromVec(V2, x = (Rand() % dimM));
370  V3.CopyRowFromMat(M, x);
371  AssertEqual(V2, V3);
372  }
373 
374  {
375  M.CopyColsFromVec(V);
376  V10.CopyColsFromMat(M);
377  AssertEqual(V, V10);
378  }
379 
380  {
381  M.CopyDiagFromVec(Vs);
382  Vs10.CopyDiagFromMat(M);
383  AssertEqual(Vs, Vs10);
384  }
385 
386  }
387 }
388 
389 template<typename Real> static void UnitTestRemoveRow() {
390 
391  // this is for matrix
392  for (MatrixIndexT p = 0;p< 10;p++) {
393  MatrixIndexT dimM = 10+Rand()%10, dimN = 10+Rand()%10;
394  Matrix<Real> M(dimM, dimN);
395  M.SetRandn();
396  MatrixIndexT i = Rand() % dimM; // Row to remove.
397  Matrix<Real> N(M);
398  N.RemoveRow(i);
399  for (MatrixIndexT j = 0;j < i;j++) {
400  for (MatrixIndexT k = 0;k < dimN;k++) {
401  KALDI_ASSERT(M(j, k) == N(j, k));
402  }
403  }
404  for (MatrixIndexT j = i+1;j < dimM;j++) {
405  for (MatrixIndexT k = 0;k < dimN;k++) {
406  KALDI_ASSERT(M(j, k) == N(j-1, k));
407  }
408  }
409  }
410 
411  // this is for vector
412  for (MatrixIndexT p = 0;p< 10;p++) {
413  MatrixIndexT dimM = 10+Rand()%10;
414  Vector<Real> V(dimM);
415  V.SetRandn();
416  MatrixIndexT i = Rand() % dimM; // Element to remove.
417  Vector<Real> N(V);
418  N.RemoveElement(i);
419  for (MatrixIndexT j = 0;j < i;j++) {
420  KALDI_ASSERT(V(j) == N(j));
421  }
422  for (MatrixIndexT j = i+1;j < dimM;j++) {
423  KALDI_ASSERT(V(j) == N(j-1));
424  }
425  }
426 
427 }
428 
429 
430 template<typename Real> static void UnitTestSubvector() {
431 
432  Vector<Real> V(100);
433  V.SetRandn();
434  Vector<Real> V2(100);
435  for (MatrixIndexT i = 0;i < 10;i++) {
436  SubVector<Real> tmp(V, i*10, 10);
437  SubVector<Real> tmp2(V2, i*10, 10);
438  tmp2.CopyFromVec(tmp);
439  }
440  AssertEqual(V, V2);
441 }
442 
443 // just need this for testing function below. Compute n!!
445  if (i <= 0) { return 1; } else { return i * DoubleFactorial(i - 2); }
446 }
447 
448 template <typename Real>
449 static void UnitTestSetRandn() {
450  for (MatrixIndexT i = 0; i < 5; i++) {
451  MatrixIndexT rows = 100 + Rand() % 50, cols = 100 + Rand() % 50;
452  Matrix<Real> M(rows, cols);
453  M.SetRandn();
454 
455  for (MatrixIndexT pow = 1; pow < 5; pow++) {
456  // test moments 1 through 4 of
457  // the distribution.
458  Matrix<Real> Mpow(M);
459  Mpow.ApplyPow(pow);
460  Real observed_moment = Mpow.Sum() / (rows * cols);
461  // see http://en.wikipedia.org/wiki/Normal_distribution#Moments,
462  // note that mu = 0 and sigma = 1.
463  Real expected_moment = (pow % 2 == 1 ? 0 : DoubleFactorial(pow - 1));
464  Real k = 10.0; // This is just a constant we use to give us some wiggle
465  // room before rejecting the distribution... e.g. 10 sigma,
466  // quite approximately.
467  Real allowed_deviation = k * pow / sqrt(static_cast<Real>(rows * cols));
468  // give it a bit more wiggle room for higher powers.. this is quite
469  // unscientific, it would be better to involve the absolute moments or
470  // something like that, and use one of those statistical inequalities,
471  // but it involves the gamma function and it's too much hassle to implement.
472  Real lower_bound = expected_moment - allowed_deviation,
473  upper_bound = expected_moment + allowed_deviation;
474  KALDI_ASSERT(observed_moment >= lower_bound && observed_moment <= upper_bound);
475  }
476  }
477 }
478 
479 
480 template <typename Real>
481 static void UnitTestSetRandUniform() {
482  for (MatrixIndexT i = 0; i < 5; i++) {
483  MatrixIndexT rows = 200 + Rand() % 50, cols = 200 + Rand() % 50;
484  Matrix<Real> M(rows, cols);
485  M.SetRandUniform();
486 
487  M.Add(-0.5); // we'll be testing the central moments, so
488  // center it around zero first.
489  // Got these moments from http://mathworld.wolfram.com/UniformDistribution.html
490  Vector<Real> central_moments(5);
491  central_moments(0) = 0.0;
492  central_moments(1) = 0.0;
493  central_moments(2) = 1.0 / 12; // times (b - a)^2, which equals 1.
494  central_moments(3) = 0.0;
495  central_moments(4) = 1.0 / 80; // times (b - a)^4, which equals 1.
496 
497  for (MatrixIndexT pow = 1; pow < central_moments.Dim(); pow++) {
498  Matrix<Real> Mpow(M);
499  Mpow.ApplyPow(pow);
500  Real observed_moment = Mpow.Sum() / (rows * cols);
501  // see http://en.wikipedia.org/wiki/Normal_distribution#Moments,
502  // note that mu = 0 and sigma = 1.
503  Real expected_moment = central_moments(pow);
504  Real k = 10.0; // This is just a constant we use to give us some wiggle
505  // room before rejecting the distribution... e.g. 10 sigma,
506  // quite approximately.
507  Real allowed_deviation = k / sqrt(static_cast<Real>(rows * cols));
508  Real lower_bound = expected_moment - allowed_deviation,
509  upper_bound = expected_moment + allowed_deviation;
510  KALDI_ASSERT(observed_moment >= lower_bound && observed_moment <= upper_bound);
511  }
512  }
513 }
514 
515 
516 template <typename Real>
517 static void UnitTestSimpleForVec() { // testing some simple operators on vector
518 
519  for (MatrixIndexT i = 0; i < 5; i++) {
520  Vector<Real> V(100), V1(100), V2(100), V3(100), V4(100);
521  V.SetRandn();
522  if (i % 2 == 0) {
523  V1.SetZero();
524  V1.Add(1.0);
525  } else {
526  V1.Set(1.0);
527  }
528 
529  V2.CopyFromVec(V);
530  V3.CopyFromVec(V1);
531  V2.InvertElements();
532  V3.DivElements(V);
533  V4.CopyFromVec(V3);
534  V4.AddVecDivVec(1.0, V1, V, 0.0);
535  AssertEqual(V3, V2);
536  AssertEqual(V4, V3);
537  V4.MulElements(V);
538  AssertEqual(V4, V1);
539  V3.AddVecVec(1.0, V, V2, 0.0);
540  AssertEqual(V3, V1);
541 
542  Vector<Real> V5(V), V6(V1), V7(V1), V8(V);
543  V5.AddVec(1.0, V);
544  V8.Scale(2.0);
545  V6.AddVec2(1.0, V);
546  V7.AddVecVec(1.0, V, V, 1.0);
547  AssertEqual(V6, V7);
548  AssertEqual(V5, V8);
549  }
550 
551  for (MatrixIndexT i = 0; i < 5; i++) {
552  std::vector<MatrixIndexT> sizes;
553  sizes.push_back(16);
554  sizes.push_back(128);
555  for(int i = 0; i < sizes.size(); i++) {
556  MatrixIndexT dimM = sizes[i] + Rand() % 10, dimN = sizes[i] + Rand() % 10;
557  Matrix<Real> M(dimM, dimN);
558  M.SetRandn();
559  Vector<Real> Vr(dimN), Vc(dimM);
560  Vr.AddRowSumMat(0.4, M);
561  Vr.AddRowSumMat(0.3, M, 0.5); // note: 0.3 + 0.4*0.5 = 0.5.
562  Vc.AddColSumMat(0.4, M);
563  Vc.AddColSumMat(0.3, M, 0.5); // note: 0.3 + 0.4*0.5 = 0.5.
564  Vr.Scale(2.0);
565  Vc.Scale(2.0);
566  KALDI_LOG << Vr;
567  KALDI_LOG << Vc;
568 
569  Vector<Real> V2r(dimN), V2c(dimM);
570  for (MatrixIndexT k = 0; k < dimM; k++) {
571  V2r.CopyRowFromMat(M, k);
572  KALDI_ASSERT(fabs(V2r.Sum() - Vc(k)) < 0.01);
573  }
574  for (MatrixIndexT j = 0; j < dimN; j++) {
575  V2c.CopyColFromMat(M, j);
576  KALDI_ASSERT(fabs(V2c.Sum() - Vr(j)) < 0.01);
577  }
578  }
579  }
580 
581  for (MatrixIndexT i = 0; i < 5; i++) {
582  Vector<Real> V(100), V1(100), V2(100);
583  V.SetRandn();
584 
585  V1.CopyFromVec(V);
586  V1.ApplyExp();
587  Real a = V.LogSumExp();
588  V2.Set(Exp(V.LogSumExp()));
589  V1.DivElements(V2);
590  V2.CopyFromVec(V);
591 
592  Real b = V.ApplySoftMax();
593  AssertEqual(V1, V);
594  AssertEqual(a, b);
595 
596  V.ApplyLog();
597  Real c = V2.ApplyLogSoftMax();
598  AssertEqual(V2, V);
599  AssertEqual(a, c);
600  }
601 
602  for (MatrixIndexT i = 0; i < 5; i++) {
603  MatrixIndexT dimV = 10 + Rand() % 10;
604  Real p = 0.5 + RandUniform() * 4.5;
605  Vector<Real> V(dimV), V1(dimV), V2(dimV);
606  V.SetRandn();
607  V1.AddVecVec(1.0, V, V, 0.0); // V1:=V.*V.
608  V2.CopyFromVec(V1);
609  AssertEqual(V1.Norm(p), V2.Norm(p));
610  AssertEqual(sqrt(V1.Sum()), V.Norm(2.0));
611  }
612 
613  for (MatrixIndexT i = 0; i < 5; i++) {
614  MatrixIndexT dimV = 10 + Rand() % 10;
615  Real p = RandUniform() * 1.0e-5;
616  Vector<Real> V(dimV);
617  V.Set(p);
618  KALDI_ASSERT(V.IsZero(p));
619  KALDI_ASSERT(!V.IsZero(p*0.9));
620  }
621 
622  // Test ApplySoftMax() for matrix.
623  Matrix<Real> M(10, 10);
624  M.SetRandn();
625  M.ApplySoftMax();
626  KALDI_ASSERT( fabs(1.0 - M.Sum()) < 0.01);
627 
628 }
629 
630 template<typename Real>
631 static void UnitTestVectorMax() {
632  int32 dimM = 1 + Rand() % 10;
633  Vector<Real> V(dimM);
634  V.SetRandn();
635  Real m = V(0);
636  for (int32 i = 1; i < dimM; i++) m = std::max(m, V(i));
637  KALDI_ASSERT(m == V.Max());
638  MatrixIndexT i;
639  KALDI_ASSERT(m == V.Max(&i));
640  KALDI_ASSERT(m == V(i));
641 }
642 
643 template<typename Real>
644 static void UnitTestVectorMin() {
645  int32 dimM = 1 + Rand() % 10;
646  Vector<Real> V(dimM);
647  V.SetRandn();
648  Real m = V(0);
649  for (int32 i = 1; i < dimM; i++) m = std::min(m, V(i));
650  KALDI_ASSERT(m == V.Min());
651  MatrixIndexT i;
652  KALDI_ASSERT(m == V.Min(&i));
653  KALDI_ASSERT(m == V(i));
654 }
655 
656 template<typename Real>
657 static void UnitTestReplaceValue(){
658  // for vector
659  MatrixIndexT dim = 10 + Rand() % 2;
660  Real orig = 0.1 * (Rand() % 100), changed = 0.1 * (Rand() % 50);
661  Vector<Real> V(dim);
662  V.SetRandn();
663  V(dim / 2) = orig;
664  Vector<Real> V1(V);
665  for (MatrixIndexT i = 0; i < V1.Dim(); i ++) {
666  if (V1(i) == orig) V1(i) = changed;
667  }
668  V.ReplaceValue(orig, changed);
669  AssertEqual(V, V1);
670 }
671 
672 
673 template<typename Real>
674 static void UnitTestNorm() { // test some simple norm properties: scaling. also ApproxEqual test.
675 
676  for (MatrixIndexT p = 0; p < 10; p++) {
677  Real scalar = RandGauss();
678  if (scalar == 0.0) continue;
679  if (scalar < 0) scalar *= -1.0;
680  MatrixIndexT dimM = 10 + Rand() % 10, dimN = 10 + Rand() % 10;
681  Matrix<Real> M(dimM, dimN);
682  M.SetRandn();
683  SpMatrix<Real> S(dimM);
684  S.SetRandn();
685  Vector<Real> V(dimN);
686  V.SetRandn();
687 
688  Real Mnorm = M.FrobeniusNorm(),
689  Snorm = S.FrobeniusNorm(),
690  Vnorm1 = V.Norm(1.0),
691  Vnorm2 = V.Norm(2.0),
692  Vnorm3 = V.Norm(3.0);
693  M.Scale(scalar);
694  S.Scale(scalar);
695  V.Scale(scalar);
696  KALDI_ASSERT(ApproxEqual(M.FrobeniusNorm(), Mnorm*scalar));
697  KALDI_ASSERT(ApproxEqual(S.FrobeniusNorm(), Snorm*scalar));
698  KALDI_ASSERT(ApproxEqual(V.Norm(1.0), Vnorm1 * scalar));
699  KALDI_ASSERT(ApproxEqual(V.Norm(2.0), Vnorm2 * scalar));
700  KALDI_ASSERT(ApproxEqual(V.Norm(3.0), Vnorm3 * scalar));
701 
705  SpMatrix<Real> S2(S); S2.Scale(1.1); KALDI_ASSERT(!S.ApproxEqual(S2)); KALDI_ASSERT(S.ApproxEqual(S2, 0.15));
706  Matrix<Real> M2(M); M2.Scale(1.1); KALDI_ASSERT(!M.ApproxEqual(M2)); KALDI_ASSERT(M.ApproxEqual(M2, 0.15));
707  Vector<Real> V2(V); V2.Scale(1.1); KALDI_ASSERT(!V.ApproxEqual(V2)); KALDI_ASSERT(V.ApproxEqual(V2, 0.15));
708  }
709 }
710 
711 
712 template<typename Real>
713 static void UnitTestCopyRows() {
714  for (MatrixIndexT p = 0; p < 10; p++) {
715  MatrixIndexT num_rows1 = 10 + Rand() % 10,
716  num_rows2 = 10 + Rand() % 10,
717  num_cols = 10 + Rand() % 10;
718  Matrix<Real> M(num_rows1, num_cols);
719  M.SetRandn();
720 
721  Matrix<Real> N1(num_rows2, num_cols),
722  N2(num_rows2, num_cols), O(num_rows2, num_cols);
723  std::vector<int32> reorder(num_rows2);
724  std::vector<const Real*> reorder_src(num_rows2,
725  static_cast<const Real*>(NULL));
726  for (int32 i = 0; i < num_rows2; i++) {
727  reorder[i] = -1 + (Rand() % (num_rows1 + 1));
728  if (reorder[i] != -1)
729  reorder_src[i] = M.RowData(reorder[i]);
730  }
731 
732  N1.CopyRows(M, &(reorder[0]));
733  N2.CopyRows(&(reorder_src[0]));
734 
735  for (int32 i = 0; i < num_rows2; i++)
736  for (int32 j = 0; j < num_cols; j++)
737  if (reorder[i] < 0) O(i, j) = 0;
738  else O(i, j) = M(reorder[i], j);
739 
740  AssertEqual(N1, O);
741  AssertEqual(N2, O);
742  }
743 }
744 
745 template<typename Real>
746 static void UnitTestCopyToRows() {
747  for (MatrixIndexT p = 0; p < 10; p++) {
748  MatrixIndexT num_rows1 = 10 + Rand() % 10,
749  num_rows2 = 10 + Rand() % 10,
750  num_cols = 10 + Rand() % 10;
751  Matrix<Real> M(num_rows1, num_cols);
752  M.SetRandn();
753 
754  Matrix<Real> N(num_rows2, num_cols), O(num_rows2, num_cols);
755  std::vector<Real*> reorder_dst(num_rows1,
756  static_cast<Real*>(NULL));
757  unordered_map<MatrixIndexT, bool> used_index;
758  for (int32 i = 0; i < num_rows1; i++) {
759  MatrixIndexT index = -1 + (Rand() % (num_rows2 + 1));
760  if (used_index.find(index) == used_index.end()) {
761  used_index[index] = true;
762  } else {
763  index = -1;
764  }
765  if (index != -1) {
766  reorder_dst[i] = N.RowData(index);
767  for (int32 j = 0; j < num_cols; j++)
768  O(index, j) = M(i, j);
769  }
770  }
771 
772  M.CopyToRows(&(reorder_dst[0]));
773 
774  AssertEqual(N, O);
775  }
776 }
777 
778 template<typename Real>
779 static void UnitTestAddRows() {
780  for (MatrixIndexT p = 0; p < 10; p++) {
781  MatrixIndexT num_rows1 = 10 + Rand() % 10,
782  num_rows2 = 10 + Rand() % 10,
783  num_cols = 10 + Rand() % 10;
784  Matrix<Real> M(num_rows1, num_cols);
785  M.SetRandn();
786 
787  Matrix<Real> N1(num_rows2, num_cols),
788  N2(num_rows2, num_cols), O(num_rows2, num_cols);
789  std::vector<int32> reorder(num_rows2);
790  std::vector<const Real*> reorder_src(num_rows2,
791  static_cast<const Real*>(NULL));
792  for (int32 i = 0; i < num_rows2; i++) {
793  reorder[i] = -1 + (Rand() % (num_rows1 + 1));
794  if (reorder[i] != -1)
795  reorder_src[i] = M.RowData(reorder[i]);
796  }
797 
798  Real alpha =
799  static_cast<Real>((Rand() % num_rows2)) / static_cast<Real>(num_rows1);
800 
801  N1.AddRows(alpha, M, &(reorder[0]));
802  N2.AddRows(alpha, &(reorder_src[0]));
803 
804  for (int32 i = 0; i < num_rows2; i++) {
805  if (reorder[i] != -1) {
806  for (int32 j = 0; j < num_cols; j++) {
807  O(i, j) += alpha * M(reorder[i], j);
808  }
809  }
810  }
811 
812  AssertEqual(N1, O);
813  AssertEqual(N2, O);
814  }
815 }
816 
817 template<typename Real>
818 static void UnitTestAddToRows() {
819  for (MatrixIndexT p = 0; p < 10; p++) {
820  MatrixIndexT num_rows1 = 10 + Rand() % 10,
821  num_rows2 = 10 + Rand() % 10,
822  num_cols = 10 + Rand() % 10;
823  Matrix<Real> M(num_rows1, num_cols);
824  M.SetRandn();
825 
826  Real alpha =
827  static_cast<Real>((Rand() % num_rows2)) / static_cast<Real>(num_rows1);
828 
829  Matrix<Real> N(num_rows2, num_cols), O(num_rows2, num_cols);
830  std::vector<Real*> reorder_dst(num_rows1, static_cast<Real*>(NULL));
831  unordered_map<MatrixIndexT, bool> used_index;
832  for (int32 i = 0; i < num_rows1; i++) {
833  MatrixIndexT index = -1 + (Rand() % (num_rows2 + 1));
834  if (used_index.find(index) == used_index.end()) {
835  used_index[index] = true;
836  } else {
837  index = -1;
838  }
839  if (index != -1) {
840  reorder_dst[i] = N.RowData(index);
841  for (int32 j = 0; j < num_cols; j++)
842  O(index, j) += alpha * M(i, j);
843  }
844  }
845 
846  M.AddToRows(alpha, &(reorder_dst[0]));
847 
848  AssertEqual(N, O);
849  }
850 }
851 
852 template<typename Real>
853 static void UnitTestCopyCols() {
854  for (MatrixIndexT p = 0; p < 10; p++) {
855  MatrixIndexT num_cols1 = 10 + Rand() % 10,
856  num_cols2 = 10 + Rand() % 10,
857  num_rows = 10 + Rand() % 10;
858  Matrix<Real> M(num_rows, num_cols1);
859  M.SetRandn();
860 
861  Matrix<Real> N(num_rows, num_cols2), O(num_rows, num_cols2);
862  std::vector<int32> reorder(num_cols2);
863  for (int32 i = 0; i < num_cols2; i++)
864  reorder[i] = -1 + (Rand() % (num_cols1 + 1));
865 
866  N.CopyCols(M, &(reorder[0]));
867 
868  for (int32 i = 0; i < num_rows; i++)
869  for (int32 j = 0; j < num_cols2; j++)
870  if (reorder[j] < 0) O(i, j) = 0;
871  else O(i, j) = M(i, reorder[j]);
872  AssertEqual(N, O);
873  }
874 }
875 
876 
877 template<typename Real>
878 static void UnitTestSimpleForMat() { // test some simple operates on all kinds of matrix
879 
880  for (MatrixIndexT p = 0; p < 10; p++) {
881  // for FrobeniousNorm() function
882  MatrixIndexT dimM = 10 + Rand() % 10, dimN = 10 + Rand() % 10;
883  Matrix<Real> M(dimM, dimN);
884  M.SetRandn();
885  {
886  Matrix<Real> N(M);
887  Real a = M.LogSumExp(), b = N.ApplySoftMax();
888  AssertEqual(a, b);
889  AssertEqual(1.0, N.Sum());
890  }
891  {
892  Matrix<Real> N(M);
893  N.Add(2.0);
894  for (MatrixIndexT m = 0; m < dimM; m++)
895  for (MatrixIndexT n = 0; n < dimN; n++)
896  N(m, n) -= 2.0;
897  AssertEqual(M, N);
898  }
899 
900  Matrix<Real> N(M), M1(M);
901  M1.MulElements(M);
902  Real tmp1 = sqrt(M1.Sum());
903  Real tmp2 = N.FrobeniusNorm();
904  KALDI_ASSERT(std::abs(tmp1 - tmp2) < 0.00001);
905 
906  // for LargestAbsElem() function
907  Vector<Real> V(dimM);
908  for (MatrixIndexT i = 0; i < dimM; i++) {
909  for (MatrixIndexT j = 0; j < dimN; j++) {
910  M(i, j) = std::abs(M(i, j));
911  }
912  std::sort(M.RowData(i), M.RowData(i) + dimN);
913  V(i) = M(i, dimN - 1);
914  }
915  std::sort(V.Data(), V.Data() + dimM);
916  KALDI_ASSERT(std::abs(V(dimM - 1) - N.LargestAbsElem()) < 0.00001);
917  }
918 
919  SpMatrix<Real> x(3);
920  x.SetZero();
921 
922  std::stringstream ss;
923 
924  ss << "DP 3\n";
925  ss << "4.6863" << '\n';
926  ss << "3.7062 4.6032" << '\n';
927  ss << "3.4160 3.7256 5.2474" << '\n';
928 
929  ss >> x;
930  KALDI_ASSERT(x.IsPosDef() == true); // test IsPosDef() function
931 
932  TpMatrix<Real> y(3);
933  y.SetZero();
934  y.Cholesky(x);
935 
936 
937  // test sp-matrix's LogPosDefDet() function
938  Matrix<Real> B(x);
939  Real tmp;
940  Real *DetSign = &tmp;
941  KALDI_ASSERT(std::abs(B.LogDet(DetSign) - x.LogPosDefDet()) < 0.00001);
942 
943  for (MatrixIndexT p = 0; p < 10; p++) { // test for sp and tp matrix's AddSp() and AddTp() function
944  MatrixIndexT dimM = 10 + Rand() % 10;
945  SpMatrix<Real> S(dimM), S1(dimM);
946  TpMatrix<Real> T(dimM), T1(dimM);
947  S.SetRandn();
948  S1.SetRandn();
949  T.SetRandn();
950  T1.SetRandn();
951  Matrix<Real> M(S), M1(S1), N(T), N1(T1);
952 
953  S.AddSp(1.0, S1);
954  T.AddTp(1.0, T1);
955  M.AddMat(1.0, M1);
956  N.AddMat(1.0, N1);
957  Matrix<Real> S2(S);
958  Matrix<Real> T2(T);
959 
960  AssertEqual(S2, M);
961  AssertEqual(T2, N);
962  }
963 
964  for (MatrixIndexT i = 0; i < 10; i++) { // test for sp matrix's AddVec2() function
965  MatrixIndexT dimM = 10 + Rand() % 10;
966  SpMatrix<Real> M(dimM);
967  Vector<Real> V(dimM);
968 
970  SpMatrix<double> Md(M);
971  V.SetRandn();
972  SpMatrix<Real> Sorig(M);
973  M.AddVec2(0.5, V);
974  Md.AddVec2(static_cast<Real>(0.5), V);
975  for (MatrixIndexT i = 0; i < dimM; i++)
976  for (MatrixIndexT j = 0; j < dimM; j++) {
977  KALDI_ASSERT(std::abs(M(i, j) - (Sorig(i, j)+0.5*V(i)*V(j))) < 0.001);
978  KALDI_ASSERT(std::abs(Md(i, j) - (Sorig(i, j)+0.5*V(i)*V(j))) < 0.001);
979  }
980  }
981 }
982 
983 
984 template<typename Real> static void UnitTestRow() {
985 
986  for (MatrixIndexT p = 0;p< 10;p++) {
987  MatrixIndexT dimM = 10+Rand()%10, dimN = 10+Rand()%10;
988  Matrix<Real> M(dimM, dimN);
990 
991  MatrixIndexT i = Rand() % dimM; // Row to get.
992 
993  Vector<Real> V(dimN);
994  V.CopyRowFromMat(M, i); // get row.
995  for (MatrixIndexT k = 0;k < dimN;k++) {
996  AssertEqual(M(i, k), V(k));
997  }
998 
999  {
1000  SpMatrix<Real> S(dimN);
1001  InitRandNonsingular(&S);
1002  Vector<Real> v1(dimN), v2(dimN);
1003  Matrix<Real> M(S);
1004  MatrixIndexT dim2 = Rand() % dimN;
1005  v1.CopyRowFromSp(S, dim2);
1006  v2.CopyRowFromMat(M, dim2);
1007  AssertEqual(v1, v2);
1008  }
1009 
1010  MatrixIndexT j = Rand() % dimN; // Col to get.
1011  Vector<Real> W(dimM);
1012  W.CopyColFromMat(M, j); // get row.
1013  for (MatrixIndexT k = 0;k < dimM;k++) {
1014  AssertEqual(M(k, j), W(k));
1015  }
1016 
1017  }
1018 }
1019 
1020 template<typename Real> static void UnitTestAxpy() {
1021 
1022  for (MatrixIndexT i = 0;i< 10;i++) {
1023  MatrixIndexT dimM = 10+Rand()%10, dimN = 10+Rand()%10;
1024  Matrix<Real> M(dimM, dimN), N(dimM, dimN), O(dimN, dimM);
1025 
1027  Matrix<Real> Morig(M);
1028  M.AddMat(0.5, N);
1029  for (MatrixIndexT i = 0;i < dimM;i++)
1030  for (MatrixIndexT j = 0;j < dimN;j++)
1031  KALDI_ASSERT(std::abs(M(i, j) - (Morig(i, j)+0.5*N(i, j))) < 0.1);
1032  M.CopyFromMat(Morig);
1033  M.AddMat(0.5, O, kTrans);
1034  for (MatrixIndexT i = 0;i < dimM;i++)
1035  for (MatrixIndexT j = 0;j < dimN;j++)
1036  KALDI_ASSERT(std::abs(M(i, j) - (Morig(i, j)+0.5*O(j, i))) < 0.1);
1037  {
1038  float f = 0.5 * (float) (Rand() % 3);
1039  Matrix<Real> N(dimM, dimM);
1040  InitRandNonsingular(&N);
1041 
1042  Matrix<Real> N2(N);
1043  Matrix<Real> N3(N);
1044  N2.AddMat(f, N2, kTrans);
1045  N3.AddMat(f, N, kTrans);
1046  AssertEqual(N2, N3); // check works same with self as arg.
1047  }
1048  }
1049 }
1050 
1051 template<typename Real> static void UnitTestCopySp() {
1052  // Checking that the various versions of copying
1053  // matrix to SpMatrix work the same in the symmetric case.
1054  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1055  int32 dim = 5 + Rand() % 10;
1056  SpMatrix<Real> S(dim), T(dim);
1057  S.SetRandn();
1058  Matrix<Real> M(S);
1060  AssertEqual(S, T);
1061  T.SetZero();
1062  T.CopyFromMat(M, kTakeMean);
1063  AssertEqual(S, T);
1064  T.SetZero();
1065  T.CopyFromMat(M, kTakeLower);
1066  AssertEqual(S, T);
1067  T.SetZero();
1068  T.CopyFromMat(M, kTakeUpper);
1069  AssertEqual(S, T);
1070  }
1071 }
1072 
1073 
1074 template<typename Real> static void UnitTestPower() {
1075  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1076  // this is for matrix-pow
1077  MatrixIndexT dimM = 10 + Rand() % 10;
1078  Matrix<Real> M(dimM, dimM), N(dimM, dimM);
1079  M.SetRandn();
1080  N.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); // N:=M*M^T.
1081  SpMatrix<Real> S(dimM);
1082  S.CopyFromMat(N); // symmetric so should not crash.
1083  S.ApplyPow(0.5);
1084  S.ApplyPow(2.0);
1085  M.CopyFromSp(S);
1086  AssertEqual(M, N);
1087 
1088  // this is for vector-pow
1089  MatrixIndexT dimV = 10 + Rand() % 10;
1090  Vector<Real> V(dimV), V1(dimV), V2(dimV);
1091  V.SetRandn();
1092  V1.AddVecVec(1.0, V, V, 0.0); // V1:=V.*V.
1093  V2.CopyFromVec(V1);
1094  V2.ApplyPow(0.5);
1095  V2.ApplyPow(2.0);
1096  AssertEqual(V1, V2);
1097  }
1098 }
1099 
1100 template<typename Real> static void UnitTestPowerAbs() {
1101  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1102  MatrixIndexT dimV = 10 + Rand() % 10;
1103  Vector<Real> V(dimV), V1(dimV), V2(dimV);
1104  V.SetRandn();
1105  V1.AddVecVec(1.0, V, V, 0.0); // V1:=V.*V.
1106  V2.CopyFromVec(V1);
1107  KALDI_LOG << V1;
1108  V2.ApplyPowAbs(0.5);
1109  KALDI_LOG << V2;
1110  V2.ApplyPowAbs(2.0);
1111  KALDI_LOG << V2;
1112  AssertEqual(V1, V2);
1113  }
1114 }
1115 
1116 
1117 template<typename Real> static void UnitTestHeaviside() {
1118  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1119  MatrixIndexT dimM = 10 + Rand() % 10, dimN = 10 + Rand() % 10;
1120  Matrix<Real> M(dimM, dimN), N(dimM, dimN);
1121  M.SetRandn();
1122  N = M;
1123  N.ApplyHeaviside();
1124  for (MatrixIndexT r = 0; r < dimM; r++) {
1125  for (MatrixIndexT c = 0; c < dimN; c++) {
1126  Real x = M(r, c), y = N(r, c);
1127  if (x < 0.0) KALDI_ASSERT(y == 0.0);
1128  if (x > 0.0) KALDI_ASSERT(y == 1.0);
1129  if (x == 0.0) { KALDI_ASSERT(y >= 0.0 && y <= 1.0); }
1130  }
1131  }
1132  }
1133 }
1134 
1135 
1136 template<typename Real> static void UnitTestAddOuterProductPlusMinus() {
1137  for (MatrixIndexT iter = 0; iter < 10; iter++) {
1138  MatrixIndexT dimM = 10 + Rand() % 10;
1139  MatrixIndexT dimN = 10 + Rand() % 10;
1140  Matrix<Real> M(dimM, dimN), Plus(dimM, dimN), Minus(dimM, dimN),
1141  M2(dimM, dimN);
1142  Vector<Real> v1(dimM), v2(dimN);
1143 
1144  for (MatrixIndexT i = 0; i < 5; i++) {
1145  v1.SetRandn();
1146  v2.SetRandn();
1147  Real alpha = 0.333 * ((Rand() % 10) - 5);
1148  M.AddVecVec(alpha, v1, v2);
1149 
1150  AddOuterProductPlusMinus(alpha, v1, v2, &Plus, &Minus);
1151  M2.SetZero();
1152  M2.AddMat(-1.0, Minus);
1153  M2.AddMat(1.0, Plus);
1154  AssertEqual(M, M2);
1155  KALDI_ASSERT(Minus.Min() >= 0);
1156  KALDI_ASSERT(Plus.Min() >= 0);
1157  }
1158  }
1159 }
1160 
1161 template<typename Real> static void UnitTestSger() {
1162  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1163  MatrixIndexT dimM = 10 + Rand() % 10;
1164  MatrixIndexT dimN = 10 + Rand() % 10;
1165  Matrix<Real> M(dimM, dimN), M2(dimM, dimN);
1166  Vector<Real> v1(dimM); v1.SetRandn();
1167  Vector<Real> v2(dimN); v2.SetRandn();
1168  Vector<double> v1d(v1), v2d(v2);
1169  M.AddVecVec(1.0f, v1, v2);
1170  M2.AddVecVec(1.0f, v1, v2);
1171  for (MatrixIndexT m = 0;m < dimM;m++)
1172  for (MatrixIndexT n = 0;n < dimN;n++) {
1173  KALDI_ASSERT(M(m, n) - v1(m)*v2(n) < 0.01);
1174  KALDI_ASSERT(M(m, n) - M2(m, n) < 0.01);
1175  }
1176  }
1177 }
1178 
1179 
1180 
1181 template<typename Real> static void UnitTestDeterminant() { // also tests matrix axpy and IsZero() and TraceOfProduct{, T}
1182  for (MatrixIndexT iter = 0;iter < 5;iter++) { // First test the 2 det routines are the same
1183  int dimM = 10 + Rand() % 10;
1184  Matrix<Real> M(dimM, dimM), N(dimM, dimM);
1185  InitRandNonsingular(&M);
1186  N.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); // N:=M*M^T.
1187  for (MatrixIndexT i = 0;i < (MatrixIndexT)dimM;i++) N(i, i) += 0.0001; // Make sure numerically +ve det-- can by chance be almost singular the way we initialized it (I think)
1188  SpMatrix<Real> S(dimM);
1189  S.CopyFromMat(N); // symmetric so should not crash.
1190  Real logdet = S.LogPosDefDet();
1191  Real logdet2, logdet3, sign2, sign3;
1192  logdet2 = N.LogDet(&sign2);
1193  logdet3 = S.LogDet(&sign3);
1194  KALDI_ASSERT(sign2 == 1.0 && sign3 == 1.0 && std::abs(logdet2-logdet) < 0.1 && std::abs(logdet2 - logdet3) < 0.1);
1195  Matrix<Real> tmp(dimM, dimM); tmp.SetZero();
1196  tmp.AddMat(1.0, N);
1197  tmp.AddMat(-1.0, N, kTrans);
1198  // symmetric so tmp should be zero.
1199  if ( ! tmp.IsZero(1.0e-04)) {
1200  printf("KALDI_ERR: matrix is not zero\n");
1201  KALDI_LOG << tmp;
1202  KALDI_ASSERT(0);
1203  }
1204 
1205  Real a = TraceSpSp(S, S), b = TraceMatMat(N, N), c = TraceMatMat(N, N, kTrans);
1206  KALDI_ASSERT(std::abs(a-b) < 0.1 && std::abs(b-c) < 0.1);
1207  }
1208 }
1209 
1210 
1211 template<typename Real> static void UnitTestDeterminantSign() {
1212 
1213  for (MatrixIndexT iter = 0;iter < 20;iter++) { // First test the 2 det routines are the same
1214  int dimM = 10 + Rand() % 10;
1215  Matrix<Real> M(dimM, dimM), N(dimM, dimM);
1216  InitRandNonsingular(&M);
1217  N.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); // N:=M*M^T.
1218  for (MatrixIndexT i = 0;i < (MatrixIndexT)dimM;i++) N(i, i) += 0.0001; // Make sure numerically +ve det-- can by chance be almost singular the way we initialized it (I think)
1219  SpMatrix<Real> S(dimM);
1220  S.CopyFromMat(N); // symmetric so should not crash.
1221  Real logdet = S.LogPosDefDet();
1222  Real logdet2, logdet3, sign2, sign3;
1223  logdet2 = N.LogDet(&sign2);
1224  logdet3 = S.LogDet(&sign3);
1225  KALDI_ASSERT(sign2 == 1.0 && sign3 == 1.0 && std::abs(logdet2-logdet) < 0.01 && std::abs(logdet2 - logdet3) < 0.01);
1226 
1227  MatrixIndexT num_sign_changes = Rand() % 5;
1228  for (MatrixIndexT change = 0; change < num_sign_changes; change++) {
1229  // Change sign of S's det by flipping one eigenvalue, and N by flipping one row.
1230  {
1231  Matrix<Real> M(S);
1232  Matrix<Real> U(dimM, dimM), Vt(dimM, dimM);
1233  Vector<Real> s(dimM);
1234  M.Svd(&s, &U, &Vt); // SVD: M = U diag(s) Vt
1235  s(Rand() % dimM) *= -1;
1236  U.MulColsVec(s);
1237  M.AddMatMat(1.0, U, kNoTrans, Vt, kNoTrans, 0.0);
1238  S.CopyFromMat(M);
1239  }
1240  // change sign of N:
1241  N.Row(Rand() % dimM).Scale(-1.0);
1242  }
1243 
1244  // add in a scaling factor too.
1245  Real tmp = 1.0 + ((Rand() % 5) * 0.01);
1246  Real logdet_factor = dimM * Log(tmp);
1247  N.Scale(tmp);
1248  S.Scale(tmp);
1249 
1250  Real logdet4, logdet5, sign4, sign5;
1251  logdet4 = N.LogDet(&sign4);
1252  logdet5 = S.LogDet(&sign5);
1253  AssertEqual(logdet4, logdet+logdet_factor, 0.01);
1254  AssertEqual(logdet5, logdet+logdet_factor, 0.01);
1255  if (num_sign_changes % 2 == 0) {
1256  KALDI_ASSERT(sign4 == 1);
1257  KALDI_ASSERT(sign5 == 1);
1258  } else {
1259  KALDI_ASSERT(sign4 == -1);
1260  KALDI_ASSERT(sign5 == -1);
1261  }
1262  }
1263 }
1264 
1265 template<typename Real> static void UnitTestSpVec() {
1266  // Test conversion back and forth between SpMatrix and Vector.
1267  for (MatrixIndexT iter = 0;iter < 1;iter++) {
1268  MatrixIndexT dimM =10; // 20 + Rand()%10;
1269  SpMatrix<Real> A(dimM), B(dimM);
1270  SubVector<Real> vec(A);
1271  B.CopyFromVec(vec);
1272  AssertEqual(A, B);
1273  }
1274 }
1275 
1276 
1277 template<typename Real> static void UnitTestTraceProduct() {
1278  for (MatrixIndexT iter = 0;iter < 5;iter++) { // First test the 2 det routines are the same
1279  int dimM = 10 + Rand() % 10, dimN = 10 + Rand() % 10;
1280  Matrix<Real> M(dimM, dimN), N(dimM, dimN);
1281 
1282  M.SetRandn();
1283  N.SetRandn();
1284  Matrix<Real> Nt(N, kTrans);
1285  Real a = TraceMatMat(M, Nt), b = TraceMatMat(M, N, kTrans);
1286  printf("m = %d, n = %d\n", dimM, dimN);
1287  KALDI_LOG << a << " " << b;
1288  KALDI_ASSERT(std::abs(a-b) < 0.1);
1289  }
1290 }
1291 
1292 template<typename Real> static void UnitTestSvd() {
1293  MatrixIndexT Base = 3, Rand_ = 2, Iter = 25;
1294  for (MatrixIndexT iter = 0;iter < Iter;iter++) {
1295  MatrixIndexT dimM = Base + Rand() % Rand_, dimN = Base + Rand() % Rand_;
1296  Matrix<Real> M(dimM, dimN);
1297  Matrix<Real> U(dimM, std::min(dimM, dimN)), Vt(std::min(dimM, dimN), dimN);
1298  Vector<Real> s(std::min(dimM, dimN));
1299  M.SetRandn();
1300  if (iter < 2) KALDI_LOG << "M " << M;
1301  Matrix<Real> M2(dimM, dimN); M2.CopyFromMat(M);
1302  M.Svd(&s, &U, &Vt);
1303  if (iter < 2) {
1304  KALDI_LOG << " s " << s;
1305  KALDI_LOG << " U " << U;
1306  KALDI_LOG << " Vt " << Vt;
1307  }
1308  MatrixIndexT min_dim = std::min(dimM, dimN);
1309  Matrix<Real> S(min_dim, min_dim);
1310  S.CopyDiagFromVec(s);
1311  Matrix<Real> Mtmp(dimM, dimN);
1312  Mtmp.SetZero();
1313  Mtmp.AddMatMatMat(1.0, U, kNoTrans, S, kNoTrans, Vt, kNoTrans, 0.0);
1314  AssertEqual(Mtmp, M2);
1315  }
1316 }
1317 
1318 template<typename Real> static void UnitTestSvdBad() {
1319  MatrixIndexT N = 20;
1320  {
1321  Matrix<Real> M(N, N);
1322  // M.Set(1591.3614306764898);
1323  M.Set(1.0);
1324  M(0, 0) *= 1.000001;
1325  Matrix<Real> U(N, N), V(N, N);
1326  Vector<Real> l(N);
1327  M.Svd(&l, &U, &V);
1328  }
1329  SpMatrix<Real> M(N);
1330  for(MatrixIndexT i =0; i < N; i++)
1331  for(MatrixIndexT j = 0; j <= i; j++)
1332  M(i, j) = 1591.3614306764898;
1333  M(0, 0) *= 1.00001;
1334  M(10, 10) *= 1.00001;
1335  Matrix<Real> U(N, N);
1336  Vector<Real> l(N);
1337  M.SymPosSemiDefEig(&l, &U);
1338 }
1339 
1340 
1341 template<typename Real> static void UnitTestSvdZero() {
1342  MatrixIndexT Base = 3, Rand_ = 2, Iter = 30;
1343  for (MatrixIndexT iter = 0;iter < Iter;iter++) {
1344  MatrixIndexT dimM = Base + Rand() % Rand_, dimN = Base + Rand() % Rand_; // M>=N.
1345  Matrix<Real> M(dimM, dimN);
1346  Matrix<Real> U(dimM, dimM), Vt(dimN, dimN); Vector<Real> v(std::min(dimM, dimN));
1347  if (iter%2 == 0) M.SetZero();
1348  else M.Unit();
1349  if (iter < 2) KALDI_LOG << "M " << M;
1350  Matrix<Real> M2(dimM, dimN); M2.CopyFromMat(M);
1351  bool ans = M.Svd(&v, &U, &Vt);
1352  KALDI_ASSERT(ans); // make sure works with zero matrix.
1353  }
1354 }
1355 
1356 
1357 
1358 
1359 
1360 template<typename Real> static void UnitTestSvdNodestroy() {
1361  MatrixIndexT Base = 3, Rand_ = 2, Iter = 25;
1362  for (MatrixIndexT iter = 0;iter < Iter;iter++) {
1363  MatrixIndexT dimN = Base + Rand() % Rand_, dimM = dimN + Rand() % Rand_; // M>=N, as required by JAMA Svd.
1364  MatrixIndexT minsz = std::min(dimM, dimN);
1365  Matrix<Real> M(dimM, dimN);
1366  Matrix<Real> U(dimM, minsz), Vt(minsz, dimN); Vector<Real> v(minsz);
1367  M.SetRandn();
1368  if (iter < 2) KALDI_LOG << "M " << M;
1369  M.Svd(&v, &U, &Vt);
1370  if (iter < 2) {
1371  KALDI_LOG << " v " << v;
1372  KALDI_LOG << " U " << U;
1373  KALDI_LOG << " Vt " << Vt;
1374  }
1375 
1376  for (MatrixIndexT it = 0;it < 2;it++) {
1377  Matrix<Real> Mtmp(minsz, minsz);
1378  for (MatrixIndexT i = 0;i < v.Dim();i++) Mtmp(i, i) = v(i);
1379  Matrix<Real> Mtmp2(minsz, dimN);
1380  Mtmp2.AddMatMat(1.0, Mtmp, kNoTrans, Vt, kNoTrans, 0.0);
1381  Matrix<Real> Mtmp3(dimM, dimN);
1382  Mtmp3.AddMatMat(1.0, U, kNoTrans, Mtmp2, kNoTrans, 0.0);
1383  for (MatrixIndexT i = 0;i < Mtmp.NumRows();i++) {
1384  for (MatrixIndexT j = 0;j < Mtmp.NumCols();j++) {
1385  KALDI_ASSERT(std::abs(Mtmp3(i, j) - M(i, j)) < 0.0001);
1386  }
1387  }
1388 
1389  SortSvd(&v, &U, &Vt); // and re-check...
1390  for (MatrixIndexT i = 0; i + 1 < v.Dim(); i++) // check SortSvd is working.
1391  KALDI_ASSERT(std::abs(v(i+1)) <= std::abs(v(i)));
1392  }
1393  }
1394 }
1395 
1396 
1397 
1398 template<typename Real> static void UnitTestSvdJustvec() { // Making sure gives same answer if we get just the vector, not the eigs.
1399  MatrixIndexT Base = 10, Rand_ = 5, Iter = 25;
1400  for (MatrixIndexT iter = 0;iter < Iter;iter++) {
1401  MatrixIndexT dimM = Base + Rand() % Rand_, dimN = Base + Rand() % Rand_; // M>=N.
1402  MatrixIndexT minsz = std::min(dimM, dimN);
1403 
1404  Matrix<Real> M(dimM, dimN);
1405  Matrix<Real> U(dimM, minsz), Vt(minsz, dimN); Vector<Real> v(minsz);
1406  M.Svd(&v, &U, &Vt);
1407  Vector<Real> v2(minsz);
1408  M.Svd(&v2);
1409  AssertEqual(v, v2);
1410  }
1411 }
1412 
1413 template<typename Real> static void UnitTestEigSymmetric() {
1414 
1415  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1416  MatrixIndexT dimM = 20 + Rand()%10;
1417  SpMatrix<Real> S(dimM);
1418  S.SetRandn();
1419  Matrix<Real> M(S); // copy to regular matrix.
1420  Matrix<Real> P(dimM, dimM);
1421  Vector<Real> real_eigs(dimM), imag_eigs(dimM);
1422  M.Eig(&P, &real_eigs, &imag_eigs);
1423  KALDI_ASSERT(imag_eigs.Sum() == 0.0);
1424  // should have M = P D P^T
1425  Matrix<Real> tmp(P); tmp.MulColsVec(real_eigs); // P * eigs
1426  Matrix<Real> M2(dimM, dimM);
1427  M2.AddMatMat(1.0, tmp, kNoTrans, P, kTrans, 0.0); // M2 = tmp * Pinv = P * eigs * P^T
1428  AssertEqual(M, M2); // check reconstruction worked.
1429  }
1430 }
1431 
1432 template<typename Real> static void UnitTestEig() {
1433 
1434  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1435  MatrixIndexT dimM = 1 + iter;
1436  /* if (iter < 10)
1437  dimM = 1 + Rand() % 6;
1438  else
1439  dimM = 5 + Rand()%10; */
1440  Matrix<Real> M(dimM, dimM);
1441  InitRandNonsingular(&M);
1442  Matrix<Real> P(dimM, dimM);
1443  Vector<Real> real_eigs(dimM), imag_eigs(dimM);
1444  M.Eig(&P, &real_eigs, &imag_eigs);
1445 
1446  { // Check that the eigenvalues match up with the determinant.
1447  BaseFloat logdet_check = 0.0, logdet = M.LogDet();
1448  for (MatrixIndexT i = 0; i < dimM ; i++)
1449  logdet_check += 0.5 * Log(real_eigs(i)*real_eigs(i) + imag_eigs(i)*imag_eigs(i));
1450  AssertEqual(logdet_check, logdet);
1451  }
1452  Matrix<Real> Pinv(P);
1453  Pinv.Invert();
1454  Matrix<Real> D(dimM, dimM);
1455  CreateEigenvalueMatrix(real_eigs, imag_eigs, &D);
1456 
1457  // check that M = P D P^{-1}.
1458  Matrix<Real> tmp(dimM, dimM);
1459  tmp.AddMatMat(1.0, P, kNoTrans, D, kNoTrans, 0.0); // tmp = P * D
1460  Matrix<Real> M2(dimM, dimM);
1461  M2.AddMatMat(1.0, tmp, kNoTrans, Pinv, kNoTrans, 0.0); // M2 = tmp * Pinv = P * D * Pinv.
1462 
1463  { // print out some stuff..
1464  Matrix<Real> MP(dimM, dimM);
1465  MP.AddMatMat(1.0, M, kNoTrans, P, kNoTrans, 0.0);
1466  Matrix<Real> PD(dimM, dimM);
1467  PD.AddMatMat(1.0, P, kNoTrans, D, kNoTrans, 0.0);
1468 
1469  Matrix<Real> PinvMP(dimM, dimM);
1470  PinvMP.AddMatMat(1.0, Pinv, kNoTrans, MP, kNoTrans, 0.0);
1471  AssertEqual(MP, PD);
1472  }
1473 
1474  AssertEqual(M, M2); // check reconstruction worked.
1475  }
1476 }
1477 
1478 
1479 template<typename Real> static void UnitTestEigSp() {
1480  // Test the Eig function with SpMatrix.
1481  // We make sure to test pathological cases, that have
1482  // either large zero eigenspaces, or two large
1483  // eigenspaces with the same absolute value but +ve
1484  // and -ve. Also zero matrix.
1485 
1486  for (MatrixIndexT iter = 0; iter < 100; iter++) {
1487  MatrixIndexT dimM = 1 + (Rand() % 10);
1488  SpMatrix<Real> S(dimM);
1489 
1490  switch (iter % 5) {
1491  case 0: // zero matrix.
1492  break;
1493  case 1: // general random symmetric matrix.
1494  InitRandNonsingular(&S);
1495  break;
1496  default:
1497  { // start with a random matrix; do decomposition; change the eigenvalues to
1498  // try to cover the problematic cases; reconstruct.
1499  InitRandNonsingular(&S);
1500  Vector<Real> s(dimM); Matrix<Real> P(dimM, dimM);
1501  S.Eig(&s, &P);
1502  // We on purpose create a problematic case where
1503  // some eigs are either zero or share a value (+ve or -ve)
1504  // with some other eigenvalue.
1505  for (MatrixIndexT i = 0; i < dimM; i++) {
1506  if (Rand() % 10 == 0) s(i) = 0; // set that eig to zero.
1507  else if (Rand() % 10 < 2) {
1508  // set that eig to some other randomly chosen eig,
1509  // times random sign.
1510  s(i) = (Rand()%2 == 0 ? 1 : -1) * s(Rand() % dimM);
1511  }
1512  }
1513  // Reconstruct s from the eigenvalues "made problematic."
1514  S.AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
1515  Real *data = s.Data();
1516  std::sort(data, data+dimM);
1517  KALDI_LOG << "Real eigs are: " << s;
1518 
1519  }
1520  }
1521  Vector<Real> s(dimM); Matrix<Real> P(dimM, dimM);
1522  S.Eig(&s, &P);
1523  KALDI_LOG << "Found eigs are: " << s;
1524  SpMatrix<Real> S2(dimM);
1525  S2.AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
1526  {
1527  SpMatrix<Real> diff(S);
1528  diff.AddSp(-1.0, S2);
1529  Vector<Real> s(dimM); Matrix<Real> P(dimM, dimM);
1530  diff.Eig(&s, &P);
1531  KALDI_LOG << "Eigs of difference are " << s;
1532  }
1533  KALDI_ASSERT(S.ApproxEqual(S2, 1.0e-03f));
1534  }
1535 }
1536 
1537 template <typename Real>
1539  SpMatrix<Real> S(transM == kTrans ? M.NumCols() : M.NumRows());
1540  S.SetUnit();
1541  S.AddMat2(-1.0, M, transM, 1.0);
1542  Real max = 0.0;
1543  for (MatrixIndexT i = 0; i < S.NumRows(); i++)
1544  for (MatrixIndexT j = 0; j <= i; j++)
1545  max = std::max(max, std::abs(S(i, j)));
1546  return max;
1547 }
1548 
1549 template<typename Real>
1550 static Real NonDiagonalness(const SpMatrix<Real> &S) {
1551  Real max_diag = 0.0, max_offdiag = 0.0;
1552  for (MatrixIndexT i = 0; i < S.NumRows(); i++)
1553  for (MatrixIndexT j = 0; j <= i; j++) {
1554  if (i == j) { max_diag = std::max(max_diag, std::abs(S(i, j))); }
1555  else { max_offdiag = std::max(max_offdiag, std::abs(S(i, j))); }
1556  }
1557  if (max_diag == 0.0) {
1558  if (max_offdiag == 0.0) return 0.0; // perfectly diagonal.
1559  else return 1.0; // perfectly non-diagonal.
1560  } else {
1561  return max_offdiag / max_diag;
1562  }
1563 }
1564 
1565 
1566 template<typename Real>
1567 static Real NonUnitness(const SpMatrix<Real> &S) {
1568  SpMatrix<Real> tmp(S.NumRows());
1569  tmp.SetUnit();
1570  tmp.AddSp(-1.0, S);
1571  Real max = 0.0;
1572  for (MatrixIndexT i = 0; i < tmp.NumRows(); i++)
1573  for (MatrixIndexT j = 0; j <= i; j++)
1574  max = std::max(max, std::abs(tmp(i, j)));
1575  return max;
1576 }
1577 
1578 template<typename Real>
1579 static void UnitTestTridiagonalize() {
1580 
1581  {
1582  float tmp[5];
1583  tmp[4] = 1.0;
1584  cblas_Xspmv(1, 0.0, tmp+2,
1585  tmp+1, 1, 0.0, tmp+4, 1);
1586  KALDI_ASSERT(tmp[4] == 0.0);
1587  }
1588  for (MatrixIndexT i = 0; i < 4; i++) {
1589  MatrixIndexT dim = 40 + Rand() % 4;
1590  // We happened to find out that a 16x16 matrix of 27's causes problems for
1591  // Tridiagonalize.
1592  if (i == 0 || i == 1)
1593  dim = 16;
1594  SpMatrix<Real> S(dim), S2(dim), R(dim), S3(dim);
1595  Matrix<Real> Q(dim, dim);
1596  InitRandNonsingular(&S);
1597  // Very small or large scaling is challenging to qr due to squares that
1598  // could go out of range.
1599  if (Rand() % 3 == 0)
1600  S.Scale(1.0e-15);
1601  else if (Rand() % 2 == 0)
1602  S.Scale(1.0e+15);
1603  if (i == 0 || i == 1) {
1604  Matrix<Real> temp(dim, dim);
1605  if (i == 0)
1606  temp.Set(27.0);
1607  else
1608  temp.Set(-1.61558713e-27);
1609  S.CopyFromMat(temp);
1610  }
1611  SpMatrix<Real> T(S);
1612  T.Tridiagonalize(&Q);
1613  KALDI_LOG << "S trace " << S.Trace() << ", T trace " << T.Trace();
1614  // KALDI_LOG << S << "\n" << T;
1615  AssertEqual(S.Trace(), T.Trace());
1616  // Also test Trace().
1617  Real ans = 0.0;
1618  for (MatrixIndexT j = 0; j < dim; j++) ans += T(j, j);
1619  AssertEqual(ans, T.Trace());
1620  if (S.LogDet() > -50.0) {
1621  // don't check logdet equality if original logdet is very negative- could
1622  // be singular.
1623  AssertEqual(T.LogDet(), S.LogDet());
1624  }
1625  R.AddMat2(1.0, Q, kNoTrans, 0.0);
1626  KALDI_LOG << "Non-unit-ness of R is " << NonUnitness(R);
1627  KALDI_ASSERT(R.IsUnit(0.01)); // Check Q is orthogonal.
1628  S2.AddMat2Sp(1.0, Q, kTrans, T, 0.0);
1629  S3.AddMat2Sp(1.0, Q, kNoTrans, S, 0.0);
1630  //KALDI_LOG << "T is " << T;
1631  //KALDI_LOG << "S is " << S;
1632  //KALDI_LOG << "S2 (should be like S) is " << S2;
1633  //KALDI_LOG << "S3 (should be like T) is " << S3;
1634  AssertEqual(S, S2);
1635  AssertEqual(T, S3);
1636  }
1637 }
1638 
1639 template<typename Real>
1641 
1642  {
1643  float tmp[5];
1644  tmp[4] = 1.0;
1645  // cblas_sspmv(CblasRowMajor, CblasLower, 1, 0.0, tmp+2,
1646  // tmp+1, 1, 0.0, tmp+4, 1);
1647  cblas_Xspmv(1, 0.0, tmp+2,
1648  tmp+1, 1, 0.0, tmp+4, 1);
1649 
1650  KALDI_ASSERT(tmp[4] == 0.0);
1651  }
1652  for (MatrixIndexT i = 0; i < 4; i++) {
1653  MatrixIndexT dim = 50 + Rand() % 4;
1654  SpMatrix<Real> S(dim), S2(dim), R(dim), S3(dim), S4(dim);
1655  Matrix<Real> Q(dim, dim);
1656  InitRandNonsingular(&S);
1657  SpMatrix<Real> T(S);
1658  T.Tridiagonalize(&Q);
1659  KALDI_LOG << "S trace " << S.Trace() << ", T trace " << T.Trace();
1660  // KALDI_LOG << S << "\n" << T;
1661  AssertEqual(S.Trace(), T.Trace());
1662  // Also test Trace().
1663  Real ans = 0.0;
1664  for (MatrixIndexT j = 0; j < dim; j++) ans += T(j, j);
1665  AssertEqual(ans, T.Trace());
1666  AssertEqual(T.LogDet(), S.LogDet());
1667  R.AddMat2(1.0, Q, kNoTrans, 0.0);
1668  KALDI_LOG << "Non-unit-ness of R after tridiag is " << NonUnitness(R);
1669  KALDI_ASSERT(R.IsUnit(0.001)); // Check Q is orthogonal.
1670  S2.AddMat2Sp(1.0, Q, kTrans, T, 0.0);
1671  S3.AddMat2Sp(1.0, Q, kNoTrans, S, 0.0);
1672  //KALDI_LOG << "T is " << T;
1673  //KALDI_LOG << "S is " << S;
1674  //KALDI_LOG << "S2 (should be like S) is " << S2;
1675  //KALDI_LOG << "S3 (should be like T) is " << S3;
1676  AssertEqual(S, S2);
1677  AssertEqual(T, S3);
1678  SpMatrix<Real> T2(T);
1679  T2.Qr(&Q);
1680  R.AddMat2(1.0, Q, kNoTrans, 0.0);
1681  KALDI_LOG << "Non-unit-ness of R after QR is " << NonUnitness(R);
1682  KALDI_ASSERT(R.IsUnit(0.001)); // Check Q is orthogonal.
1683  AssertEqual(T.Trace(), T2.Trace());
1684  KALDI_ASSERT(T2.IsDiagonal());
1685  AssertEqual(T.LogDet(), T2.LogDet());
1686  S4.AddMat2Sp(1.0, Q, kTrans, T2, 0.0);
1687  //KALDI_LOG << "S4 (should be like S) is " << S4;
1688  AssertEqual(S, S4);
1689  }
1690 }
1691 
1692 
1693 template<typename Real> static void UnitTestMmul() {
1694  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1695  MatrixIndexT dimM = 20 + Rand()%10, dimN = 20 + Rand()%10, dimO = 20 + Rand()%10; // dims between 10 and 20.
1696  // MatrixIndexT dimM = 2, dimN = 3, dimO = 4;
1697  Matrix<Real> A(dimM, dimN), B(dimN, dimO), C(dimM, dimO);
1698  A.SetRandn();
1699  B.SetRandn();
1700  //
1701  // KALDI_LOG <<"a = " << A;
1702  // KALDI_LOG<<"B = " << B;
1703  C.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0); // C = A * B.
1704  // KALDI_LOG << "c = " << C;
1705  for (MatrixIndexT i = 0;i < dimM;i++) {
1706  for (MatrixIndexT j = 0;j < dimO;j++) {
1707  double sum = 0.0;
1708  for (MatrixIndexT k = 0;k < dimN;k++) {
1709  sum += A(i, k) * B(k, j);
1710  }
1711  KALDI_ASSERT(std::abs(sum - C(i, j)) < 0.0001);
1712  }
1713  }
1714  }
1715 }
1716 
1717 
1718 template<typename Real> static void UnitTestMmulSym() {
1719 
1720  // Test matrix multiplication on symmetric matrices.
1721  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1722  MatrixIndexT dimM = 20 + Rand()%10;
1723 
1724  Matrix<Real> A(dimM, dimM), B(dimM, dimM), C(dimM, dimM), tmp(dimM, dimM), tmp2(dimM, dimM);
1725  SpMatrix<Real> sA(dimM), sB(dimM), sC(dimM), stmp(dimM);
1726 
1727  A.SetRandn();
1728  B.SetRandn();
1729  C.SetRandn();
1730  // Make A, B, C symmetric.
1731  tmp.CopyFromMat(A); A.AddMat(1.0, tmp, kTrans);
1732  tmp.CopyFromMat(B); B.AddMat(1.0, tmp, kTrans);
1733  tmp.CopyFromMat(C); C.AddMat(1.0, tmp, kTrans);
1734 
1735  sA.CopyFromMat(A);
1736  sB.CopyFromMat(B);
1737  sC.CopyFromMat(C);
1738 
1739 
1740  tmp.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0); // tmp = A * B.
1741  tmp2.AddSpSp(1.0, sA, sB, 0.0); // tmp = sA*sB.
1742  AssertEqual(tmp, tmp2);
1743  tmp2.AddSpSp(1.0, sA, sB, 0.0); // tmp = sA*sB.
1744  AssertEqual(tmp, tmp2);
1745  }
1746 }
1747 
1748 
1749 template<typename Real> static void UnitTestAddVecVec() {
1750  for (int32 i = 0; i < 20; i++) {
1751  int32 dimM = 5 + Rand() % 10, dimN = 5 + Rand() % 10;
1752 
1753  Matrix<Real> M(dimM, dimN);
1754  M.SetRandn();
1755  Matrix<Real> N(M);
1756  Vector<float> v(dimM), w(dimN);
1757  v.SetRandn();
1758  w.SetRandn();
1759  float alpha = 0.2 * (Rand() % 10);
1760  M.AddVecVec(alpha, v, w);
1761  for (int32 j = 0; j < 20; j++) {
1762  int32 dimX = Rand() % dimM, dimY = Rand() % dimN;
1763  AssertEqual(M(dimX, dimY),
1764  N(dimX, dimY) + alpha * v(dimX) * w(dimY));
1765  }
1766  }
1767 }
1768 
1769 
1770 template<typename Real> static void UnitTestVecmul() {
1771  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1772  MatrixTransposeType trans = (iter % 2 == 0 ? kTrans : kNoTrans);
1773  MatrixIndexT dimM = 20 + Rand()%10, dimN = 20 + Rand()%10; // dims between 10 and 20.
1774  Real alpha = 0.333, beta = 0.5;
1775  Matrix<Real> A(dimM, dimN);
1776  if (trans == kTrans) A.Transpose();
1777  A.SetRandn();
1778  Vector<Real> x(dimM), y(dimN);
1779  //x.SetRandn();
1780  y.SetRandn();
1781  Vector<Real> orig_x(x), x2(x);
1782 
1783  x.AddMatVec(alpha, A, trans, y, beta); // x = A * y + beta*x.
1784  x2.AddMatSvec(alpha, A, trans, y, beta); // x = A * y + beta*x
1785 
1786  for (MatrixIndexT i = 0; i < dimM; i++) {
1787  double sum = beta * orig_x(i);
1788  for (MatrixIndexT j = 0; j < dimN; j++) {
1789  if (trans == kNoTrans) {
1790  sum += alpha * A(i, j) * y(j);
1791  } else {
1792  sum += alpha * A(j, i) * y(j);
1793  }
1794  }
1795  KALDI_ASSERT(std::abs(sum - x(i)) < 0.0001);
1796  }
1797  AssertEqual(x, x2);
1798  }
1799 
1800 }
1801 
1802 template<typename Real> static void UnitTestInverse() {
1803  for (MatrixIndexT iter = 0;iter < 10;iter++) {
1804  MatrixIndexT dimM = 20 + Rand()%10;
1805  Matrix<Real> A(dimM, dimM), B(dimM, dimM), C(dimM, dimM);
1806  InitRandNonsingular(&A);
1807  B.CopyFromMat(A);
1808  B.Invert();
1809 
1810  C.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0); // C = A * B.
1811 
1812 
1813  for (MatrixIndexT i = 0;i < dimM;i++)
1814  for (MatrixIndexT j = 0;j < dimM;j++)
1815  KALDI_ASSERT(std::abs(C(i, j) - (i == j?1.0:0.0)) < 0.1);
1816  }
1817 }
1818 
1819 
1820 
1821 
1822 template<typename Real> static void UnitTestMulElements() {
1823  for (MatrixIndexT iter = 0; iter < 5; iter++) {
1824  MatrixIndexT dimM = 20 + Rand()%10, dimN = 20 + Rand()%10;
1825  Matrix<Real> A(dimM, dimN), B(dimM, dimN), C(dimM, dimN);
1826  A.SetRandn();
1827  B.SetRandn();
1828 
1829  C.CopyFromMat(A);
1830  C.MulElements(B); // C = A .* B (in Matlab, for example).
1831 
1832  for (MatrixIndexT i = 0;i < dimM;i++)
1833  for (MatrixIndexT j = 0;j < dimN;j++)
1834  KALDI_ASSERT(std::abs(C(i, j) - (A(i, j)*B(i, j))) < 0.0001);
1835  }
1836 }
1837 
1838 template<typename Real> static void UnitTestDotprod() {
1839  for (MatrixIndexT iter = 0;iter < 5;iter++) {
1840  MatrixIndexT dimM = 200 + Rand()%100;
1841  Vector<Real> v(dimM), w(dimM);
1842 
1843  v.SetRandn();
1844  w.SetRandn();
1845  Vector<double> wd(w);
1846 
1847  Real f = VecVec(w, v), f2 = VecVec(wd, v), f3 = VecVec(v, wd);
1848  Real sum = 0.0;
1849  for (MatrixIndexT i = 0;i < dimM;i++) sum += v(i)*w(i);
1850  KALDI_ASSERT(std::abs(f-sum) < 0.001);
1851  KALDI_ASSERT(std::abs(f2-sum) < 0.001);
1852  KALDI_ASSERT(std::abs(f3-sum) < 0.001);
1853  }
1854 }
1855 
1856 
1857 template<class Real>
1859  int32 num_rows = mat->NumRows(), num_cols = mat->NumCols(),
1860  stride = mat->Stride();
1861  BaseFloat not_a_number = nan(" "); // nan is from <cmath>
1862  for (int32 r = 0; r + 1 < num_rows; r++) {
1863  for (int32 j = num_cols; j < stride; j++) {
1864  if (RandInt(0, 1) == 0)
1865  (mat->RowData(r))[j] = not_a_number;
1866  else
1867  (mat->RowData(r))[j] = RandGauss() * 1.5e+31;
1868  }
1869  }
1870 }
1871 
1872 
1878 template<typename Real>
1880  for (size_t i = 0; i < 10; i++) {
1881  MatrixIndexT num_rows = Rand() % 10, num_cols = Rand() % 10;
1882  if (num_rows * num_cols == 0) num_rows = num_cols = 0;
1883  MatrixStrideType src_stride_type = (Rand() % 2 == 0) ?
1885  // Always pick the other stride type.
1886  MatrixStrideType resize_stride_type = (src_stride_type == kDefaultStride) ?
1888  Matrix<Real> src(num_rows, num_cols, kSetZero, src_stride_type);
1889  Matrix<Real> clone(src);
1890  PlaceNansInGaps(&clone);
1891  src.Resize(num_rows, num_cols, kCopyData, resize_stride_type);
1892  PlaceNansInGaps(&src);
1893  AssertEqual(src, clone);
1894  }
1895 }
1896 
1897 
1898 template<typename Real>
1899 static void UnitTestResize() {
1900  for (size_t i = 0; i < 10; i++) {
1901  MatrixIndexT dimM1 = Rand() % 10, dimN1 = Rand() % 10,
1902  dimM2 = Rand() % 10, dimN2 = Rand() % 10;
1903  if (dimM1*dimN1 == 0) dimM1 = dimN1 = 0;
1904  if (dimM2*dimN2 == 0) dimM2 = dimN2 = 0;
1905  for (MatrixIndexT j = 0; j < 3; j++) {
1906  MatrixResizeType resize_type = static_cast<MatrixResizeType>(j);
1907  Matrix<Real> M(dimM1, dimN1);
1908  M.SetRandn();
1909  Matrix<Real> Mcopy(M);
1910  Vector<Real> v(dimM1);
1911  v.SetRandn();
1912  Vector<Real> vcopy(v);
1913  SpMatrix<Real> S(dimM1);
1914  S.SetRandn();
1915  SpMatrix<Real> Scopy(S);
1916  M.Resize(dimM2, dimN2, resize_type);
1917  v.Resize(dimM2, resize_type);
1918  S.Resize(dimM2, resize_type);
1919  if (resize_type == kSetZero) {
1920  KALDI_ASSERT(S.IsZero());
1921  KALDI_ASSERT(v.Sum() == 0.0);
1922  KALDI_ASSERT(M.IsZero());
1923  } else if (resize_type == kCopyData) {
1924  for (MatrixIndexT i = 0; i < dimM2; i++) {
1925  if (i < dimM1) AssertEqual(v(i), vcopy(i));
1926  else KALDI_ASSERT(v(i) == 0);
1927  for (MatrixIndexT j = 0; j < dimN2; j++) {
1928  if (i < dimM1 && j < dimN1) AssertEqual(M(i, j), Mcopy(i, j));
1929  else AssertEqual(M(i, j), 0.0);
1930  }
1931  for (MatrixIndexT i2 = 0; i2 < dimM2; i2++) {
1932  if (i < dimM1 && i2 < dimM1) AssertEqual(S(i, i2), Scopy(i, i2));
1933  else AssertEqual(S(i, i2), 0.0);
1934  }
1935  }
1936  }
1937  }
1938  }
1939 }
1940 
1941 
1942 template<typename Real>
1943 static void UnitTestTp2Sp() {
1944  // Tests AddTp2Sp()
1945  for (MatrixIndexT iter = 0; iter < 4; iter++) {
1946  MatrixIndexT dimM = 10 + Rand()%3;
1947 
1948  TpMatrix<Real> T(dimM);
1949  T.SetRandn();
1950  SpMatrix<Real> S(dimM);
1951  S.SetRandn();
1952 
1953  Matrix<Real> M(T);
1954  for ( MatrixIndexT i = 0; i < dimM; i++)
1955  for (MatrixIndexT j = 0; j < dimM; j++) {
1956  if (j <= i) AssertEqual(T(i,j), M(i,j));
1957  else AssertEqual(M(i,j), 0.0);
1958  }
1959 
1960  SpMatrix<Real> A(dimM), B(dimM);
1961  A.AddTp2Sp(0.5, T, (iter < 2 ? kNoTrans : kTrans), S, 0.0);
1962  B.AddMat2Sp(0.5, M, (iter < 2 ? kNoTrans : kTrans), S, 0.0);
1963  AssertEqual(A, B);
1964  }
1965 }
1966 
1967 template<typename Real>
1968 static void UnitTestTp2() {
1969  // Tests AddTp2()
1970  for (MatrixIndexT iter = 0; iter < 4; iter++) {
1971  MatrixIndexT dimM = 10 + Rand()%3;
1972 
1973  TpMatrix<Real> T(dimM);
1974  T.SetRandn();
1975 
1976  Matrix<Real> M(T);
1977 
1978  SpMatrix<Real> A(dimM), B(dimM);
1979  A.AddTp2(0.5, T, (iter < 2 ? kNoTrans : kTrans), 0.0);
1980  B.AddMat2(0.5, M, (iter < 2 ? kNoTrans : kTrans), 0.0);
1981  AssertEqual(A, B);
1982  }
1983 }
1984 
1985 template<typename Real>
1986 static void UnitTestAddDiagMat2() {
1987  for (MatrixIndexT iter = 0; iter < 4; iter++) {
1988  MatrixIndexT dimM = 10 + Rand() % 3,
1989  dimN = 1 + Rand() % 4;
1990  Vector<Real> v(dimM);
1991  v.SetRandn();
1992  Vector<Real> w(v);
1993  if (iter % 2 == 1) {
1994  Matrix<Real> M(dimM, dimN);
1995  M.SetRandn();
1996  v.AddDiagMat2(0.5, M, kNoTrans, 0.3);
1997  Matrix<Real> M2(dimM, dimM);
1998  M2.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0);
1999  Vector<Real> diag(dimM);
2000  diag.CopyDiagFromMat(M2);
2001  w.Scale(0.3);
2002  w.AddVec(0.5, diag);
2003  AssertEqual(w, v);
2004  } else {
2005  Matrix<Real> M(dimN, dimM);
2006  M.SetRandn();
2007  v.AddDiagMat2(0.5, M, kTrans, 0.3);
2008  Matrix<Real> M2(dimM, dimM);
2009  M2.AddMatMat(1.0, M, kTrans, M, kNoTrans, 0.0);
2010  Vector<Real> diag(dimM);
2011  diag.CopyDiagFromMat(M2);
2012  w.Scale(0.3);
2013  w.AddVec(0.5, diag);
2014  AssertEqual(w, v);
2015  }
2016  }
2017 }
2018 
2019 template<typename Real>
2020 static void UnitTestAddDiagMatMat() {
2021  for (MatrixIndexT iter = 0; iter < 4; iter++) {
2022  BaseFloat alpha = 0.432 + Rand() % 5, beta = 0.043 + Rand() % 2;
2023  MatrixIndexT dimM = 10 + Rand() % 3,
2024  dimN = 5 + Rand() % 4;
2025  Vector<Real> v(dimM);
2026  Matrix<Real> M_orig(dimM, dimN), N_orig(dimN, dimM);
2027  M_orig.SetRandn();
2028  N_orig.SetRandn();
2029  MatrixTransposeType transM = (iter % 2 == 0 ? kNoTrans : kTrans);
2030  MatrixTransposeType transN = ((iter/2) % 2 == 0 ? kNoTrans : kTrans);
2031  Matrix<Real> M(M_orig, transM), N(N_orig, transN);
2032 
2033  v.SetRandn();
2034  Vector<Real> w(v);
2035 
2036  w.AddDiagMatMat(alpha, M, transM, N, transN, beta);
2037 
2038  {
2039  Vector<Real> w2(v);
2040  Matrix<Real> MN(dimM, dimM);
2041  MN.AddMatMat(1.0, M, transM, N, transN, 0.0);
2042  Vector<Real> d(dimM);
2043  d.CopyDiagFromMat(MN);
2044  w2.Scale(beta);
2045  w2.AddVec(alpha, d);
2046  AssertEqual(w, w2);
2047  }
2048  }
2049 }
2050 
2051 template<typename Real>
2053  for (MatrixIndexT iter = 0; iter < 100; iter++) {
2054  MatrixIndexT dimM = 4 + Rand() % 5, dimN = dimM + (Rand() % 2);
2055  Matrix<Real> M(dimM, dimN);
2056  for (MatrixIndexT i = 0; i < dimM; i++) {
2057  if (Rand() % 5 != 0) M.Row(i).SetRandn();
2058  }
2059  if (Rand() % 2 != 0) { // Multiply by a random square matrix;
2060  // keeps it low rank but will be correlated. Harder
2061  // test case.
2062  Matrix<Real> N(dimM, dimM);
2063  N.SetRandn();
2064  Matrix<Real> tmp(dimM, dimN);
2065  tmp.AddMatMat(1.0, N, kNoTrans, M, kNoTrans, 0.0);
2066  M.CopyFromMat(tmp);
2067  }
2068  M.OrthogonalizeRows();
2069  Matrix<Real> I(dimM, dimM);
2070  I.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0);
2071  KALDI_ASSERT(I.IsUnit(1.0e-05));
2072  }
2073 }
2074 
2075 template<typename Real>
2077  for (MatrixIndexT iter = 0;iter < 10;iter++) {
2078 
2079  MatrixIndexT dimA = 10 + Rand()%3;
2080  MatrixIndexT dimO = 10 + Rand()%3;
2081  Matrix<Real> Af(dimA, dimA);
2082  SpMatrix<Real> Ap(dimA);
2083  Matrix<Real> M(dimO, dimA);
2084  Matrix<Real> Of(dimO, dimO);
2085  SpMatrix<Real> Op(dimO);
2086  Matrix<Real> A_MT(dimA, dimO);
2087 
2088  for (MatrixIndexT i = 0;i < Ap.NumRows();i++) {
2089  for (MatrixIndexT j = 0; j<=i; j++) {
2090  Ap(i, j) = RandGauss();
2091  }
2092  }
2093  for (MatrixIndexT i = 0;i < M.NumRows();i++) {
2094  for (MatrixIndexT j = 0; j < M.NumCols(); j++) {
2095  M(i, j) = RandGauss();
2096  }
2097  }
2098  /*
2099  std::stringstream ss("1 2 3");
2100  ss >> Ap;
2101  ss.clear();
2102  ss.str("5 6 7 8 9 10");
2103  ss >> M;
2104  */
2105 
2106  Af.CopyFromSp(Ap);
2107  A_MT.AddMatMat(1.0, Af, kNoTrans, M, kTrans, 0.0);
2108  Of.AddMatMat(1.0, M, kNoTrans, A_MT, kNoTrans, 0.0);
2109  Op.AddMat2Sp(1.0, M, kNoTrans, Ap, 0.0);
2110 
2111 
2112  // KALDI_LOG << "A" << '\n' << Af;
2113  // KALDI_LOG << "M" << '\n' << M;
2114  // KALDI_LOG << "Op" << '\n' << Op;
2115 
2116  for (MatrixIndexT i = 0; i < dimO; i++) {
2117  for (MatrixIndexT j = 0; j<=i; j++) {
2118  KALDI_ASSERT(std::abs(Of(i, j) - Op(i, j)) < 0.0001);
2119  }
2120  }
2121 
2122  A_MT.Resize(dimO, dimA);
2123  A_MT.AddMatMat(1.0, Of, kNoTrans, M, kNoTrans, 0.0);
2124  Af.AddMatMat(1.0, M, kTrans, A_MT, kNoTrans, 1.0);
2125  Ap.AddMat2Sp(1.0, M, kTrans, Op, 1.0);
2126 
2127  // KALDI_LOG << "Ap" << '\n' << Ap;
2128  // KALDI_LOG << "Af" << '\n' << Af;
2129 
2130  for (MatrixIndexT i = 0; i < dimA; i++) {
2131  for (MatrixIndexT j = 0; j<=i; j++) {
2132  KALDI_ASSERT(std::abs(Af(i, j) - Ap(i, j)) < 0.01);
2133  }
2134  }
2135  }
2136 }
2137 
2138 
2139 template<typename Real>
2140 static void UnitTestRankNUpdate() {
2141  for (MatrixIndexT iter = 0;iter < 10;iter++) {
2142  MatrixIndexT dimA = 10 + Rand()%3;
2143  MatrixIndexT dimO = 10 + Rand()%3;
2144  Matrix<Real> Af(dimA, dimA);
2145  SpMatrix<Real> Ap(dimA);
2146  SpMatrix<Real> Ap2(dimA);
2147  Matrix<Real> M(dimO, dimA);
2148  M.SetRandn();
2149  Matrix<Real> N(M, kTrans);
2150  Af.AddMatMat(1.0, M, kTrans, M, kNoTrans, 0.0);
2151  Ap.AddMat2(1.0, M, kTrans, 0.0);
2152  Ap2.AddMat2(1.0, N, kNoTrans, 0.0);
2153  Matrix<Real> Ap_f(Ap);
2154  Matrix<Real> Ap2_f(Ap2);
2155  AssertEqual(Ap_f, Af);
2156  AssertEqual(Ap2_f, Af);
2157  }
2158 }
2159 
2160 template<typename Real> static void UnitTestSpInvert() {
2161  for (MatrixIndexT i = 0;i < 30;i++) {
2162  MatrixIndexT dimM = 6 + Rand()%20;
2163  SpMatrix<Real> M(dimM);
2164  for (MatrixIndexT i = 0;i < M.NumRows();i++)
2165  for (MatrixIndexT j = 0;j<=i;j++) M(i, j) = RandGauss();
2166  SpMatrix<Real> N(dimM);
2167  N.CopyFromSp(M);
2168  if (Rand() % 2 == 0)
2169  N.Invert();
2170  else
2171  N.InvertDouble();
2172  Matrix<Real> Mm(dimM, dimM), Nm(dimM, dimM), Om(dimM, dimM);
2173  Mm.CopyFromSp(M); Nm.CopyFromSp(N);
2174  Om.AddMatMat(1.0, Mm, kNoTrans, Nm, kNoTrans, 0.0);
2175  KALDI_ASSERT(Om.IsUnit( 0.01*dimM ));
2176  }
2177 }
2178 
2179 
2180 template<typename Real> static void UnitTestTpInvert() {
2181  for (MatrixIndexT i = 0;i < 30;i++) {
2182  MatrixIndexT dimM = 20 + Rand()%10;
2183  TpMatrix<Real> M(dimM);
2184  for (MatrixIndexT i = 0;i < M.NumRows();i++) {
2185  for (MatrixIndexT j = 0;j < i;j++) M(i, j) = RandGauss();
2186  M(i, i) = 20 * std::max((Real)0.1, (Real) RandGauss()); // make sure invertible by making it diagonally dominant (-ish)
2187  }
2188  TpMatrix<Real> N(dimM);
2189  N.CopyFromTp(M);
2190  N.Invert();
2191  TpMatrix<Real> O(dimM);
2192 
2193  Matrix<Real> Mm(dimM, dimM), Nm(dimM, dimM), Om(dimM, dimM);
2194  Mm.CopyFromTp(M); Nm.CopyFromTp(N);
2195 
2196  Om.AddMatMat(1.0, Mm, kNoTrans, Nm, kNoTrans, 0.0);
2197  KALDI_ASSERT(Om.IsUnit(0.001));
2198  }
2199 }
2200 
2201 
2202 template<typename Real> static void UnitTestLimitCondInvert() {
2203  for (MatrixIndexT i = 0;i < 10;i++) {
2204  MatrixIndexT dimM = 20 + Rand()%10;
2205  MatrixIndexT dimN = dimM + 1 + Rand()%10;
2206 
2207  SpMatrix<Real> B(dimM);
2208  Matrix<Real> X(dimM, dimN);
2209  X.SetRandn();
2210  B.AddMat2(1.0, X, kNoTrans, 0.0); // B = X*X^T -> positive definite (almost certainly), since N > M.
2211 
2212 
2213  SpMatrix<Real> B2(B);
2214  B2.LimitCond(1.0e+10, true); // Will invert.
2215 
2216  Matrix<Real> Bf(B), B2f(B2);
2217  Matrix<Real> I(dimM, dimM); I.AddMatMat(1.0, Bf, kNoTrans, B2f, kNoTrans, 0.0);
2218  KALDI_ASSERT(I.IsUnit(0.1));
2219  }
2220 }
2221 
2222 
2223 template<typename Real> static void UnitTestFloorChol() {
2224  for (MatrixIndexT i = 0;i < 10;i++) {
2225  MatrixIndexT dimM = 20 + Rand()%10;
2226 
2227 
2228  MatrixIndexT dimN = 20 + Rand()%10;
2229  Matrix<Real> X(dimM, dimN);
2230  X.SetRandn();
2231  SpMatrix<Real> B(dimM);
2232  B.AddMat2(1.0, X, kNoTrans, 0.0); // B = X*X^T -> positive semidefinite.
2233 
2234  float alpha = (Rand() % 10) + 0.5;
2235  Matrix<Real> M(dimM, dimM);
2236  M.SetRandn();
2237  SpMatrix<Real> C(dimM);
2238  C.AddMat2(1.0, M, kNoTrans, 0.0); // C:=M*M^T
2239  InitRandNonsingular(&M);
2240  C.AddMat2(1.0, M, kNoTrans, 1.0); // C+=M*M^T (after making new random M)
2241  if (i%2 == 0)
2242  C.Scale(0.001); // so it's not too much bigger than B (or it's trivial)
2243  SpMatrix<Real> BFloored(B); BFloored.ApplyFloor(C, alpha);
2244 
2245 
2246  for (MatrixIndexT j = 0;j < 10;j++) {
2247  Vector<Real> v(dimM);
2248  v.SetRandn();
2249  Real ip_b = VecSpVec(v, B, v);
2250  Real ip_a = VecSpVec(v, BFloored, v);
2251  Real ip_c = alpha * VecSpVec(v, C, v);
2252  if (i < 3) KALDI_LOG << "alpha = " << alpha << ", ip_a = " << ip_a << " ip_b = " << ip_b << " ip_c = " << ip_c <<'\n';
2253  KALDI_ASSERT(ip_a>=ip_b*0.999 && ip_a>=ip_c*0.999);
2254  }
2255  }
2256 }
2257 
2258 
2259 
2260 
2261 template<typename Real> static void UnitTestFloorUnit() {
2262  for (MatrixIndexT i = 0;i < 5;i++) {
2263  MatrixIndexT dimM = 20 + Rand()%10;
2264  MatrixIndexT dimN = 20 + Rand()%10;
2265  float floor = (Rand() % 10) - 3;
2266 
2267  Matrix<Real> M(dimM, dimN);
2268  M.SetRandn();
2269  SpMatrix<Real> B(dimM);
2270  B.AddMat2(1.0, M, kNoTrans, 0.0); // B = M*M^T -> positive semidefinite.
2271 
2272  SpMatrix<Real> BFloored(B);
2273  BFloored.ApplyFloor(floor);
2274 
2275 
2276  Vector<Real> s(dimM); Matrix<Real> P(dimM, dimM); B.SymPosSemiDefEig(&s, &P);
2277  Vector<Real> s2(dimM); Matrix<Real> P2(dimM, dimM); BFloored.SymPosSemiDefEig(&s2, &P2);
2278 
2279  KALDI_ASSERT ( (s.Min() >= floor && std::abs(s2.Min()-s.Min())<0.01) || std::abs(s2.Min()-floor)<0.01);
2280  }
2281 }
2282 
2283 
2284 template<typename Real> static void UnitTestFloorCeiling() {
2285  for (MatrixIndexT i = 0; i < 5; i++) {
2286  MatrixIndexT dimM = 10 + Rand() % 10;
2287  Vector<Real> v(dimM);
2288  v.SetRandn();
2289  Real pivot = v(5);
2290  Vector<Real> f(v), f2(v), c(v), c2(v);
2291  MatrixIndexT floored2;
2292  f.ApplyFloor(pivot, &floored2);
2293  MatrixIndexT ceiled2;
2294  c.ApplyCeiling(pivot, &ceiled2);
2295  MatrixIndexT floored = 0, ceiled = 0;
2296  for (MatrixIndexT d = 0; d < dimM; d++) {
2297  if (f2(d) < pivot) { f2(d) = pivot; floored++; }
2298  if (c2(d) > pivot) { c2(d) = pivot; ceiled++; }
2299  }
2300  AssertEqual(f, f2);
2301  AssertEqual(c, c2);
2302  KALDI_ASSERT(floored == floored2);
2303  KALDI_ASSERT(ceiled == ceiled2);
2304 
2305  // Check that the non-counting variants are equivalent to the counting
2306  // variants.
2307  Vector<Real> f3(v);
2308  f3.ApplyFloor(pivot);
2309  AssertEqual(f2, f3);
2310 
2311  Vector<Real> c3(v);
2312  c3.ApplyCeiling(pivot);
2313  AssertEqual(c2, c3);
2314  }
2315 }
2316 
2317 template<typename Real> static void UnitTestMat2Vec() {
2318  for (MatrixIndexT i = 0; i < 5; i++) {
2319  MatrixIndexT dimM = 10 + Rand() % 10;
2320 
2321  Matrix<Real> M(dimM, dimM);
2322  M.SetRandn();
2323  SpMatrix<Real> B(dimM);
2324  B.AddMat2(1.0, M, kNoTrans, 0.0); // B = M*M^T -> positive definite (since M nonsingular).
2325 
2326  Matrix<Real> P(dimM, dimM);
2327  Vector<Real> s(dimM);
2328 
2329  B.SymPosSemiDefEig(&s, &P);
2330  SpMatrix<Real> B2(dimM);
2331  B2.CopyFromSp(B);
2332  B2.Scale(0.25);
2333 
2334  // B2 <-- 2.0*B2 + 0.5 * P * diag(v) * P^T
2335  B2.AddMat2Vec(0.5, P, kNoTrans, s, 2.0); // 2.0 * 0.25 + 0.5 = 1.
2336  AssertEqual(B, B2);
2337 
2338  SpMatrix<Real> B3(dimM);
2339  Matrix<Real> PT(P, kTrans);
2340  B3.AddMat2Vec(1.0, PT, kTrans, s, 0.0);
2341  AssertEqual(B, B3);
2342  }
2343 }
2344 
2345 template<typename Real> static void UnitTestLimitCond() {
2346  for (MatrixIndexT i = 0;i < 5;i++) {
2347  MatrixIndexT dimM = 20 + Rand()%10;
2348  SpMatrix<Real> B(dimM);
2349  B(1, 1) = 10000;
2350  KALDI_ASSERT(B.LimitCond(1000) == (dimM-1));
2351  KALDI_ASSERT(std::abs(B(2, 2) - 10.0) < 0.01);
2352  KALDI_ASSERT(std::abs(B(3, 0)) < 0.001);
2353  }
2354 }
2355 
2356 template<typename Real> static void UnitTestTanh() {
2357  for (MatrixIndexT i = 0; i < 10; i++) {
2358  MatrixIndexT dimM = 5 + Rand() % 10, dimN = 5 + Rand() % 10;
2359  Matrix<Real> M(dimM, dimN), P(dimM, dimN), Q(dimM, dimN), R(dimM, dimN);
2360  M.SetRandn();
2361  P.SetRandn();
2362  Matrix<Real> N(M);
2363  for(int32 r = 0; r < dimM; r++) {
2364  for (int32 c = 0; c < dimN; c++) {
2365  Real x = N(r, c);
2366  if (x > 0.0) {
2367  x = -1.0 + 2.0 / (1.0 + Exp(-2.0 * x));
2368  } else {
2369  x = 1.0 - 2.0 / (1.0 + Exp(2.0 * x));
2370  }
2371  N(r, c) = x;
2372  Real out_diff = P(r, c), in_diff = out_diff * (1.0 - x * x);
2373  Q(r, c) = in_diff;
2374  }
2375  }
2376  M.Tanh(M);
2377  R.DiffTanh(N, P);
2378  AssertEqual(M, N);
2379  AssertEqual(Q, R);
2380  }
2381 }
2382 
2383 template<typename Real> static void UnitTestSigmoid() {
2384  for (MatrixIndexT i = 0; i < 10; i++) {
2385  MatrixIndexT dimM = 5 + Rand() % 10, dimN = 5 + Rand() % 10;
2386  Matrix<Real> M(dimM, dimN), P(dimM, dimN), Q(dimM, dimN), R(dimM, dimN);
2387  M.SetRandn();
2388  P.SetRandn();
2389  Matrix<Real> N(M);
2390  for(int32 r = 0; r < dimM; r++) {
2391  for (int32 c = 0; c < dimN; c++) {
2392  Real x = N(r, c),
2393  y = 1.0 / (1 + Exp(-x));
2394  N(r, c) = y;
2395  Real out_diff = P(r, c), in_diff = out_diff * y * (1.0 - y);
2396  Q(r, c) = in_diff;
2397  }
2398  }
2399  M.Sigmoid(M);
2400  R.DiffSigmoid(N, P);
2401  AssertEqual(M, N);
2402  AssertEqual(Q, R);
2403  }
2404 }
2405 
2406 template<typename Real> static void UnitTestSoftHinge() {
2407  for (MatrixIndexT i = 0; i < 10; i++) {
2408  MatrixIndexT dimM = 5 + Rand() % 10, dimN = 5 + Rand() % 10;
2409  Matrix<Real> M(dimM, dimN), N(dimM, dimN), O(dimM, dimN);
2410  M.SetRandn();
2411  M.Scale(20.0);
2412 
2413  for(int32 r = 0; r < dimM; r++) {
2414  for (int32 c = 0; c < dimN; c++) {
2415  Real x = M(r, c);
2416  Real &y = N(r, c);
2417  if (x > 10.0) y = x;
2418  else y = Log1p(Exp(x));
2419  }
2420  }
2421  O.SoftHinge(M);
2422  AssertEqual(N, O);
2423  }
2424 }
2425 
2426 
2427 template<typename Real> static void UnitTestSimple() {
2428  for (MatrixIndexT i = 0;i < 5;i++) {
2429  MatrixIndexT dimM = 20 + Rand()%10, dimN = 20 + Rand()%20;
2430  Matrix<Real> M(dimM, dimN);
2431  M.SetUnit();
2432  KALDI_ASSERT(M.IsUnit());
2433  KALDI_ASSERT(!M.IsZero());
2434  KALDI_ASSERT(M.IsDiagonal());
2435 
2436  SpMatrix<Real> S(dimM);
2437  S.SetRandn();
2438  Matrix<Real> N(S);
2439  KALDI_ASSERT(!N.IsDiagonal()); // technically could be diagonal, but almost infinitely unlikely.
2441  KALDI_ASSERT(!N.IsUnit());
2442  KALDI_ASSERT(!N.IsZero());
2443 
2444  M.SetZero();
2445  KALDI_ASSERT(M.IsZero());
2446  Vector<Real> V(dimM*dimN);
2447  V.SetRandn();
2448  Vector<Real> V2(V), V3(dimM*dimN);
2449  V2.ApplyExp();
2450  AssertEqual(V.Sum(), V2.SumLog());
2451  V3.ApplyLogAndCopy(V2);
2452  V2.ApplyLog();
2453  AssertEqual(V, V2);
2454  AssertEqual(V3, V2);
2455 
2456  {
2457  Vector<Real> V2(V);
2458  for (MatrixIndexT i = 0; i < V2.Dim(); i++)
2459  V2(i) = Exp(V2(i));
2460  V.ApplyExp();
2461  AssertEqual(V, V2);
2462  }
2463  {
2464  Matrix<Real> N2(N), N3(N);
2465  for (MatrixIndexT i = 0; i < N.NumRows(); i++)
2466  for (MatrixIndexT j = 0; j < N.NumCols(); j++)
2467  N2(i, j) = Exp(N2(i, j));
2468  N3.ApplyExp();
2469  AssertEqual(N2, N3);
2470  }
2471  KALDI_ASSERT(!S.IsDiagonal());
2472  KALDI_ASSERT(!S.IsUnit());
2473  N.SetUnit();
2474  S.CopyFromMat(N);
2475  KALDI_ASSERT(S.IsDiagonal());
2476  KALDI_ASSERT(S.IsUnit());
2477  N.SetZero();
2478  S.CopyFromMat(N);
2479  KALDI_ASSERT(S.IsZero());
2480  KALDI_ASSERT(S.IsDiagonal());
2481  }
2482 }
2483 
2484 
2485 
2486 template<typename Real> static void UnitTestIo() {
2487 
2488  for (MatrixIndexT i = 0;i < 5;i++) {
2489  MatrixIndexT dimM = Rand()%10 + 1;
2490  MatrixIndexT dimN = Rand()%10 + 1;
2491  bool binary = (i%2 == 0);
2492 
2493  if (i == 0) {
2494  dimM = 0;dimN = 0; // test case when both are zero.
2495  }
2496  Matrix<Real> M(dimM, dimN);
2497  M.SetRandn();
2498  Matrix<Real> N;
2499  Vector<Real> v1(dimM);
2500  v1.SetRandn();
2501  Vector<Real> v2(dimM);
2502 
2503  SpMatrix<Real> S(dimM);
2504  S.SetRandn();
2505  KALDI_LOG << "SpMatrix IO: " << S;
2506  SpMatrix<Real> T(dimM);
2507 
2508  {
2509  std::ofstream outs("tmpf", std::ios_base::out |std::ios_base::binary);
2510  InitKaldiOutputStream(outs, binary);
2511  M.Write(outs, binary);
2512  S.Write(outs, binary);
2513  v1.Write(outs, binary);
2514  M.Write(outs, binary);
2515  S.Write(outs, binary);
2516  v1.Write(outs, binary);
2517  }
2518 
2519  {
2520  bool binary_in;
2521  bool either_way = (i%2 == 0);
2522  std::ifstream ins("tmpf", std::ios_base::in | std::ios_base::binary);
2523  if (!InitKaldiInputStream(ins, &binary_in)) {
2524  KALDI_ERR << "Malformed input stream.";
2525  }
2526  N.Resize(0, 0);
2527  T.Resize(0);
2528  v2.Resize(0);
2529  N.Read(ins, binary_in, either_way);
2530  T.Read(ins, binary_in, either_way);
2531  v2.Read(ins, binary_in, either_way);
2532  if (i%2 == 0)
2533  ((MatrixBase<Real>&)N).Read(ins, binary_in, true); // add
2534  else
2535  N.Read(ins, binary_in, true);
2536  T.Read(ins, binary_in, true); // add
2537  if (i%2 == 0)
2538  ((VectorBase<Real>&)v2).Read(ins, binary_in, true); // add
2539  else
2540  v2.Read(ins, binary_in, true);
2541  }
2542  N.Scale(0.5);
2543  v2.Scale(0.5);
2544  T.Scale(0.5);
2545  AssertEqual(M, N);
2546  AssertEqual(v1, v2);
2547  AssertEqual(S, T);
2548  }
2549 }
2550 
2551 
2552 template<typename Real> static void UnitTestIoCross() { // across types.
2553 
2554  typedef typename OtherReal<Real>::Real Other; // e.g. if Real == float, Other == double.
2555  for (MatrixIndexT i = 0;i < 5;i++) {
2556  MatrixIndexT dimM = Rand()%10 + 1;
2557  MatrixIndexT dimN = Rand()%10 + 1;
2558  bool binary = (i%2 == 0);
2559  if (i == 0) {
2560  dimM = 0;dimN = 0; // test case when both are zero.
2561  }
2562  Matrix<Real> M(dimM, dimN);
2563  Matrix<Other> MO;
2564  M.SetRandn();
2565  Matrix<Real> N(dimM, dimN);
2566  Vector<Real> v(dimM);
2567  Vector<Other> vO;
2568  v.SetRandn();
2569  Vector<Real> w(dimM);
2570 
2571  SpMatrix<Real> S(dimM);
2572  SpMatrix<Other> SO;
2573  S.SetRandn();
2574  SpMatrix<Real> T(dimM);
2575 
2576  {
2577  std::ofstream outs("tmpf", std::ios_base::out |std::ios_base::binary);
2578  InitKaldiOutputStream(outs, binary);
2579 
2580  M.Write(outs, binary);
2581  S.Write(outs, binary);
2582  v.Write(outs, binary);
2583  M.Write(outs, binary);
2584  S.Write(outs, binary);
2585  v.Write(outs, binary);
2586  }
2587  {
2588  std::ifstream ins("tmpf", std::ios_base::in | std::ios_base::binary);
2589  bool binary_in;
2590  if (!InitKaldiInputStream(ins, &binary_in)) {
2591  KALDI_ERR << "Malformed input stream";
2592  }
2593 
2594  MO.Read(ins, binary_in);
2595  SO.Read(ins, binary_in);
2596  vO.Read(ins, binary_in);
2597  MO.Read(ins, binary_in, true);
2598  SO.Read(ins, binary_in, true);
2599  vO.Read(ins, binary_in, true);
2600  N.CopyFromMat(MO);
2601  T.CopyFromSp(SO);
2602  w.CopyFromVec(vO);
2603  }
2604  N.Scale(0.5);
2605  w.Scale(0.5);
2606  T.Scale(0.5);
2607  AssertEqual(M, N);
2608  AssertEqual(v, w);
2609  AssertEqual(S, T);
2610  }
2611 }
2612 
2613 
2614 template<typename Real> static void UnitTestHtkIo() {
2615 
2616  for (MatrixIndexT i = 0;i < 5;i++) {
2617  MatrixIndexT dimM = Rand()%10 + 10;
2618  MatrixIndexT dimN = Rand()%10 + 10;
2619 
2620  HtkHeader hdr;
2621  hdr.mNSamples = dimM;
2622  hdr.mSamplePeriod = 10000; // in funny HTK units-- can set it arbitrarily
2623  hdr.mSampleSize = sizeof(float)*dimN;
2624  hdr.mSampleKind = 8; // Mel spectrum.
2625 
2626  Matrix<Real> M(dimM, dimN);
2627  M.SetRandn();
2628 
2629  {
2630  std::ofstream os("tmpf", std::ios::out|std::ios::binary);
2631  WriteHtk(os, M, hdr);
2632  }
2633 
2634  Matrix<Real> N;
2635  HtkHeader hdr2;
2636  {
2637  std::ifstream is("tmpf", std::ios::in|std::ios::binary);
2638  ReadHtk(is, &N, &hdr2);
2639  }
2640 
2641  AssertEqual(M, N);
2642  KALDI_ASSERT(hdr.mNSamples == hdr2.mNSamples);
2644  KALDI_ASSERT(hdr.mSampleSize == hdr2.mSampleSize);
2645  KALDI_ASSERT(hdr.mSampleKind == hdr2.mSampleKind);
2646  }
2647 
2648  unlink("tmpf");
2649 }
2650 
2651 
2652 
2653 template<typename Real> static void UnitTestRange() { // Testing SubMatrix class.
2654 
2655  // this is for matrix-range
2656  for (MatrixIndexT i = 0;i < 5;i++) {
2657  MatrixIndexT dimM = (Rand()%10) + 10;
2658  MatrixIndexT dimN = (Rand()%10) + 10;
2659 
2660  Matrix<Real> M(dimM, dimN);
2661  M.SetRandn();
2662  MatrixIndexT dimMStart = Rand() % 5;
2663  MatrixIndexT dimNStart = Rand() % 5;
2664 
2665  MatrixIndexT dimMEnd = dimMStart + 1 + (Rand()%10); if (dimMEnd > dimM) dimMEnd = dimM;
2666  MatrixIndexT dimNEnd = dimNStart + 1 + (Rand()%10); if (dimNEnd > dimN) dimNEnd = dimN;
2667 
2668 
2669  SubMatrix<Real> sub(M, dimMStart, dimMEnd-dimMStart, dimNStart, dimNEnd-dimNStart);
2670 
2671  KALDI_ASSERT(sub.Sum() == M.Range(dimMStart, dimMEnd-dimMStart, dimNStart, dimNEnd-dimNStart).Sum());
2672 
2673  for (MatrixIndexT i = dimMStart;i < dimMEnd;i++)
2674  for (MatrixIndexT j = dimNStart;j < dimNEnd;j++)
2675  KALDI_ASSERT(M(i, j) == sub(i-dimMStart, j-dimNStart));
2676 
2677  sub.SetRandn();
2678 
2679  KALDI_ASSERT(sub.Sum() == M.Range(dimMStart, dimMEnd-dimMStart, dimNStart, dimNEnd-dimNStart).Sum());
2680 
2681  for (MatrixIndexT i = dimMStart;i < dimMEnd;i++)
2682  for (MatrixIndexT j = dimNStart;j < dimNEnd;j++)
2683  KALDI_ASSERT(M(i, j) == sub(i-dimMStart, j-dimNStart));
2684  }
2685 
2686  // this if for vector-range
2687  for (MatrixIndexT i = 0;i < 5;i++) {
2688  MatrixIndexT length = (Rand()%10) + 10;
2689 
2690  Vector<Real> V(length);
2691  V.SetRandn();
2692  MatrixIndexT lenStart = Rand() % 5;
2693 
2694  MatrixIndexT lenEnd = lenStart + 1 + (Rand()%10); if (lenEnd > length) lenEnd = length;
2695 
2696  SubVector<Real> sub(V, lenStart, lenEnd-lenStart);
2697 
2698  KALDI_ASSERT(ApproxEqual(sub.Sum(), V.Range(lenStart, lenEnd-lenStart).Sum()));
2699 
2700  for (MatrixIndexT i = lenStart;i < lenEnd;i++)
2701  KALDI_ASSERT(V(i) == sub(i-lenStart));
2702 
2703  sub.SetRandn();
2704 
2705  KALDI_ASSERT(ApproxEqual(sub.Sum(), V.Range(lenStart, lenEnd-lenStart).Sum()));
2706 
2707  for (MatrixIndexT i = lenStart;i < lenEnd;i++)
2708  KALDI_ASSERT(V(i) == sub(i-lenStart));
2709  }
2710 }
2711 
2712 template<typename Real> static void UnitTestScale() {
2713 
2714  for (MatrixIndexT i = 0;i < 5;i++) {
2715  MatrixIndexT dimM = (Rand()%10) + 10;
2716  MatrixIndexT dimN = (Rand()%10) + 10;
2717 
2718  Matrix<Real> M(dimM, dimN);
2719 
2720  Matrix<Real> N(M);
2721  float f = (float)((Rand()%10)-5);
2722  M.Scale(f);
2723  KALDI_ASSERT(M.Sum() == f * N.Sum());
2724 
2725  {
2726  // now test scale_rows
2727  M.CopyFromMat(N); // make same.
2728  Vector<Real> V(dimM);
2729  V.SetRandn();
2730  M.MulRowsVec(V);
2731  for (MatrixIndexT i = 0; i < dimM;i++)
2732  for (MatrixIndexT j = 0;j < dimN;j++)
2733  KALDI_ASSERT(M(i, j) - N(i, j)*V(i) < 0.0001);
2734  }
2735 
2736  {
2737  // now test scale_cols
2738  M.CopyFromMat(N); // make same.
2739  Vector<Real> V(dimN);
2740  V.SetRandn();
2741  M.MulColsVec(V);
2742  for (MatrixIndexT i = 0; i < dimM;i++)
2743  for (MatrixIndexT j = 0;j < dimN;j++)
2744  KALDI_ASSERT(M(i, j) - N(i, j)*V(j) < 0.0001);
2745  }
2746 
2747  }
2748 }
2749 
2750 template<typename Real> static void UnitTestMul() {
2751  for (MatrixIndexT x = 0; x<=1; x++) {
2752  MatrixTransposeType trans = (x == 1 ? kTrans: kNoTrans);
2753  for (MatrixIndexT i = 0;i < 5;i++) {
2754  float alpha = 1.0, beta =0;
2755  if (i%3 == 0) beta = 0.5;
2756  if (i%5 == 0) alpha = 0.7;
2757  MatrixIndexT dimM = (Rand()%10) + 10;
2758  Vector<Real> v(dimM);
2759  v.SetRandn();
2760  TpMatrix<Real> T(dimM);
2761  T.SetRandn();
2762  Matrix<Real> M(dimM, dimM);
2763  if (i%2 == 1)
2764  M.CopyFromTp(T, trans);
2765  else
2766  M.CopyFromTp(T, kNoTrans);
2767  Vector<Real> v2(v);
2768  Vector<Real> v3(v);
2769  v2.AddTpVec(alpha, T, trans, v, beta);
2770  if (i%2 == 1)
2771  v3.AddMatVec(alpha, M, kNoTrans, v, beta);
2772  else
2773  v3.AddMatVec(alpha, M, trans, v, beta);
2774 
2775  v.AddTpVec(alpha, T, trans, v, beta);
2776  AssertEqual(v2, v3);
2777  AssertEqual(v, v2);
2778  }
2779 
2780  for (MatrixIndexT i = 0;i < 5;i++) {
2781  float alpha = 1.0, beta =0;
2782  if (i%3 == 0) beta = 0.5;
2783  if (i%5 == 0) alpha = 0.7;
2784 
2785  MatrixIndexT dimM = (Rand()%10) + 10;
2786  Vector<Real> v(dimM);
2787  v.SetRandn();
2788  SpMatrix<Real> T(dimM);
2789  T.SetRandn();
2790  Matrix<Real> M(T);
2791  Vector<Real> v2(dimM);
2792  Vector<Real> v3(dimM);
2793  v2.AddSpVec(alpha, T, v, beta);
2794  v3.AddMatVec(alpha, M, i%2 ? kNoTrans : kTrans, v, beta);
2795  AssertEqual(v2, v3);
2796  }
2797  }
2798 }
2799 
2800 template<typename Real> static void UnitTestApplyExpSpecial() {
2801  int32 rows = RandInt(1, 10), cols = RandInt(1, 10);
2802  Matrix<Real> mat(rows, cols);
2803  mat.SetRandn();
2804  Matrix<Real> A(mat), B(mat);
2805  A.ApplyExp();
2806  B.Add(1.0);
2807  B.ApplyFloor(1.0);
2808  A.Min(B); // min of exp(x) and max(1.0, x + 1).
2809  mat.ApplyExpSpecial();
2810  KALDI_LOG << "A is: " << A;
2811  AssertEqual(mat, A);
2812 }
2813 
2814 template<typename Real> static void UnitTestInnerProd() {
2815 
2816  MatrixIndexT N = 1 + Rand() % 10;
2817  SpMatrix<Real> S(N);
2818  S.SetRandn();
2819  Vector<Real> v(N);
2820  v.SetRandn();
2821  Real prod = VecSpVec(v, S, v);
2822  Real f2=0.0;
2823  for (MatrixIndexT i = 0; i < N; i++)
2824  for (MatrixIndexT j = 0; j < N; j++) {
2825  f2 += v(i) * v(j) * S(i, j);
2826  }
2827  AssertEqual(prod, f2);
2828 }
2829 
2830 
2831 template<typename Real> static void UnitTestAddToDiag() {
2832  MatrixIndexT N = 1 + Rand() % 10;
2833  SpMatrix<Real> S(N);
2834  S.SetRandn();
2835  SpMatrix<Real> S2(S);
2836  Real x = 0.5;
2837  S.AddToDiag(x);
2838  for (MatrixIndexT i = 0; i < N; i++) S2(i, i) += x;
2839  AssertEqual(S, S2);
2840 }
2841 
2842 template<typename Real> static void UnitTestScaleDiag() {
2843 
2844  MatrixIndexT N = 1 + Rand() % 10;
2845  SpMatrix<Real> S(N);
2846  S.SetRandn();
2847  SpMatrix<Real> S2(S);
2848  S.ScaleDiag(0.5);
2849  for (MatrixIndexT i = 0; i < N; i++) S2(i, i) *= 0.5;
2850  AssertEqual(S, S2);
2851 }
2852 
2853 
2854 template<typename Real> static void UnitTestSetDiag() {
2855 
2856  MatrixIndexT N = 1 + Rand() % 10;
2857  SpMatrix<Real> S(N), T(N);
2858  S.SetUnit();
2859  S.ScaleDiag(0.5);
2860  T.SetDiag(0.5);
2861  AssertEqual(S, T);
2862 
2863 }
2864 
2865 template<typename Real> static void UnitTestTraceSpSpLower() {
2866 
2867  MatrixIndexT N = 1 + Rand() % 10;
2868  SpMatrix<Real> S(N), T(N);
2869  S.SetRandn();
2870  T.SetRandn();
2871 
2872  SpMatrix<Real> Sfast(S);
2873  Sfast.Scale(2.0);
2874  Sfast.ScaleDiag(0.5);
2875 
2876  AssertEqual(TraceSpSp(S, T), TraceSpSpLower(Sfast, T));
2877 }
2878 
2879 // also tests AddSmatMat
2880 template<typename Real> static void UnitTestAddMatSmat() {
2881  for (MatrixIndexT i = 0; i < 6; i++) {
2882  MatrixIndexT dimM = (Rand()%10) + 1,
2883  dimN = (Rand()%10 + 1),
2884  dimO = (Rand()%10 + 1);
2885  MatrixTransposeType transB = (i % 2 == 0 ? kTrans : kNoTrans),
2886  transC = (i % 3 == 0 ? kTrans : kNoTrans);
2887  Matrix<Real> A(dimM, dimN),
2888  B(dimM, dimO), C(dimO, dimN);
2889  A.SetRandn(); B.SetRandn(); C.SetRandn();
2890  if (transB == kTrans) B.Transpose();
2891  if (transC == kTrans) C.Transpose();
2892  Matrix<Real> A2(A), A3(A);
2893  BaseFloat beta = 0.333, alpha = 0.5;
2894  A.AddMatMat(alpha, B, transB, C, transC, beta);
2895  A2.AddMatSmat(alpha, B, transB, C, transC, beta);
2896  A3.AddSmatMat(alpha, B, transB, C, transC, beta);
2897  AssertEqual(A, A2);
2898  AssertEqual(A, A3);
2899  }
2900 }
2901 
2902 // Also tests AddSmat2Sp
2903 template<typename Real> static void UnitTestAddMat2Sp() {
2904  for (MatrixIndexT i = 0; i < 5; i++) {
2905  MatrixIndexT dimM = (Rand()%10) + 1,
2906  dimN = (Rand()%10 + 1);
2907  BaseFloat alpha = 0.8, beta = 0.9;
2908  SpMatrix<Real> S(dimM), T(dimN);
2909  S.SetRandn();
2910  T.SetRandn();
2911  Matrix<Real> M(dimM, dimN);
2912  M.SetRandn();
2913  MatrixTransposeType trans = (i % 2 == 1 ? kTrans: kNoTrans);
2914  if (trans == kTrans) M.Transpose();
2915  SpMatrix<Real> S2(S), S3(S);
2916  S.AddMat2Sp(alpha, M, trans, T, beta);
2917  S3.AddSmat2Sp(alpha, M, trans, T, beta);
2918 
2919  // M[trans?] * T.
2920  Matrix<Real> A(dimM, dimN);
2921  A.AddMatSp(1.0, M, trans, T, 0.0);
2922  Matrix<Real> B(dimM, dimM); // M[trans] * T * M[trans]'
2923  B.AddMatMat(1.0, A, kNoTrans, M, trans == kTrans ? kNoTrans : kTrans, 0.0);
2924  SpMatrix<Real> tmp(B);
2925  S2.Scale(beta);
2926  S2.AddSp(alpha, tmp);
2927  AssertEqual(S, S2);
2928  AssertEqual(S, S3);
2929  }
2930 }
2931 
2932 template<typename Real> static void UnitTestAddMatSelf() {
2933  MatrixIndexT dimM = (Rand() % 10) + 1;
2934  Matrix<Real> M(dimM, dimM), N(dimM, dimM);
2935  M.SetRandn();
2936  N.AddMat(1.5, M);
2937  M.AddMat(0.5, M);
2938  AssertEqual(M, N);
2939  N.AddMat(0.5, M, kTrans);
2940  M.AddMat(0.5, M, kTrans);
2941  AssertEqual(M, N);
2942 }
2943 
2944 template<typename Real> static void UnitTestAddMat2() {
2945  MatrixIndexT extra = 1;
2946  // Test AddMat2 function of SpMatrix.
2947  for (MatrixIndexT i = 0; i < 5; i++) {
2948  MatrixIndexT dimM = (Rand()%10) + extra,
2949  dimN = (Rand() % 10) + extra;
2950  Real alpha = 0.2 * (Rand() % 6),
2951  beta = 0.2 * (Rand() % 6);
2952  SpMatrix<Real> S(dimM);
2953  S.SetRandn();
2954  MatrixTransposeType trans = (i % 2 == 1 ? kTrans: kNoTrans),
2955  other_trans = (trans == kTrans ? kNoTrans : kTrans);
2956  Matrix<Real> M;
2957  if (trans == kNoTrans) M.Resize(dimM, dimN);
2958  else M.Resize(dimN, dimM);
2959  M.SetRandn();
2960 
2961  Matrix<Real> Sfull(S), Sfull2(S);
2962 
2963  S.AddMat2(alpha, M, trans, beta);
2964 
2965  Sfull.AddMatMat(alpha, M, trans, M, other_trans, beta);
2966 
2967  Sfull2.SymAddMat2(alpha, M, trans, beta);
2968 
2969  // now symmetrize.
2970  SpMatrix<Real> Sfull2_copy(Sfull2, kTakeLower);
2971  Sfull2.CopyFromSp(Sfull2_copy);
2972 
2973  Matrix<Real> Sfull3(S);
2974  AssertEqual(Sfull, Sfull3);
2975  AssertEqual(Sfull2, Sfull3);
2976  }
2977 }
2978 
2979 
2980 template<typename Real> static void UnitTestSymAddMat2() {
2981  for (int32 i = 0; i < 5; i++) {
2982  int32 dimM = 10 + Rand() % 200, dimN = 10 + Rand() % 30;
2983  KALDI_LOG << "dimM = " << dimM << ", dimN = " << dimN;
2984 
2985  Matrix<Real> M(dimM, dimM); // square matrix..
2986  Matrix<Real> N(dimM, dimN);
2987  M.SetRandn();
2988  N.SetRandn();
2989  //MatrixTransposeType trans = (i % 2 == 0 ? kTrans : kNoTrans),
2990  MatrixTransposeType trans = kTrans,
2991  other_trans = (trans == kTrans ? kNoTrans : kTrans);
2992  if (trans == kTrans) N.Transpose();
2993  KALDI_LOG << "N sum is " << N.Sum();
2994  Matrix<Real> M2(M);
2995  KALDI_LOG << "M sum is " << M.Sum();
2996 
2997  Real alpha = 0.2 * (Rand() % 6),
2998  beta = 0.2 * (Rand() % 6);
2999  //Real alpha = 0.3, beta = 1.75432;
3000  M.SymAddMat2(alpha, N, trans, beta);
3001 
3002  KALDI_LOG << "M(0, 0) is " << M(0, 0);
3003  KALDI_LOG << "M sum2 is " << M.Sum();
3004 
3005  M2.AddMatMat(alpha, N, trans, N, other_trans, beta);
3006 
3007  TpMatrix<Real> T1(M.NumRows()), T2(M2.NumRows());
3008  T1.CopyFromMat(M);
3009  T2.CopyFromMat(M2);
3010  Matrix<Real> X1(T1), X2(T2); // so we can test equality.
3011  AssertEqual(X1, X2);
3012  KALDI_ASSERT(dimM == 0 || X1.Trace() != 0 || (alpha == 0 && beta == 0));
3013  }
3014 }
3015 
3016 
3017 template<typename Real> static void UnitTestSolve() {
3018 
3019  for (MatrixIndexT i = 0;i < 5;i++) {
3020  MatrixIndexT dimM = (Rand()%10) + 10;
3021  MatrixIndexT dimN = dimM - (Rand()%3); // slightly lower-dim.
3022 
3023  SpMatrix<Real> H(dimM);
3024  Matrix<Real> M(dimM, dimN);
3025  M.SetRandn();
3026  H.AddMat2(1.0, M, kNoTrans, 0.0); // H = M M^T
3027 
3028  Vector<Real> x(dimM);
3029  x.SetRandn();
3030  Vector<Real> tmp(dimM);
3031  tmp.SetRandn();
3032  Vector<Real> g(dimM);
3033  g.AddSpVec(1.0, H, tmp, 0.0); // Limit to subspace that H is in.
3034  Vector<Real> x2(x), x3(x);
3035  SolverOptions opts2, opts3;
3036  opts2.diagonal_precondition = Rand() % 2;
3037  opts2.optimize_delta = Rand() % 2;
3038  opts3.diagonal_precondition = Rand() % 2;
3039  opts3.optimize_delta = Rand() % 2;
3040 
3041  double ans2 = SolveQuadraticProblem(H, g, opts2, &x2),
3042  ans3 = SolveQuadraticProblem(H, g, opts3, &x3);
3043 
3044  double observed_impr2 = (VecVec(x2, g) -0.5* VecSpVec(x2, H, x2)) -
3045  (VecVec(x, g) -0.5* VecSpVec(x, H, x)),
3046  observed_impr3 = (VecVec(x3, g) -0.5* VecSpVec(x3, H, x3)) -
3047  (VecVec(x, g) -0.5* VecSpVec(x, H, x));
3048  AssertEqual(observed_impr2, ans2);
3049  AssertEqual(observed_impr3, ans3);
3050  KALDI_ASSERT(ans2 >= 0);
3051  KALDI_ASSERT(ans3 >= 0);
3052  KALDI_ASSERT(std::abs(ans2 - ans3) / std::max(ans2, ans3) < 0.01);
3053  //AssertEqual(x2, x3);
3054  //AssertEqual(ans1, ans2);
3055  }
3056 
3057 
3058  for (MatrixIndexT i = 0; i < 5; i++) {
3059  MatrixIndexT dimM = (Rand() % 10) + 10;
3060  MatrixIndexT dimN = dimM - (Rand() % 3); // slightly lower-dim.
3061  MatrixIndexT dimO = (Rand() % 10) + 10;
3062 
3063  SpMatrix<Real> Q(dimM), SigmaInv(dimO);
3064  Matrix<Real> Mtmp(dimM, dimN);
3065  Mtmp.SetRandn();
3066  Q.AddMat2(1.0, Mtmp, kNoTrans, 0.0); // H = M M^T
3067 
3068  Matrix<Real> Ntmp(dimO, dimN);
3069  Ntmp.SetRandn();
3070  SigmaInv.AddMat2(1.0, Ntmp, kNoTrans, 0.0); // H = M M^T
3071 
3072  Matrix<Real> M(dimO, dimM), Y(dimO, dimM);
3073  M.SetRandn();
3074  Y.SetRandn();
3075 
3076  Matrix<Real> M2(M);
3077 
3078  SpMatrix<Real> Qinv(Q);
3079  if (Q.Cond() < 1000.0) Qinv.Invert();
3080 
3081  SolverOptions opts;
3082  opts.optimize_delta = Rand() % 2;
3083  opts.diagonal_precondition = Rand() % 2;
3084  double ans = SolveQuadraticMatrixProblem(Q, Y, SigmaInv, opts, &M2);
3085 
3086  Matrix<Real> M3(M);
3087  M3.AddMatSp(1.0, Y, kNoTrans, Qinv, 0.0);
3088  if (Q.Cond() < 1000.0) {
3089  AssertEqual(M2, M3); // This equality only holds if SigmaInv full-rank,
3090  // which is overwhelmingly likely if dimO > dimM
3091  }
3092 
3093  {
3094  Real a1 = TraceMatSpMat(M2, kTrans, SigmaInv, Y, kNoTrans),
3095  a2 = TraceMatSpMatSp(M2, kNoTrans, Q, M2, kTrans, SigmaInv),
3096  b1 = TraceMatSpMat(M, kTrans, SigmaInv, Y, kNoTrans),
3097  b2 = TraceMatSpMatSp(M, kNoTrans, Q, M, kTrans, SigmaInv),
3098  a3 = a1 - 0.5 * a2,
3099  b3 = b1 - 0.5 * b2;
3100  KALDI_ASSERT(a3 >= b3);
3101  AssertEqual(a3 - b3, ans);
3102  // KALDI_LOG << "a3 = " << a3 << ", b3 = " << b3 << ", c3 = " << c3;
3103  } // Check objf not decreased.
3104  }
3105 
3106  for (MatrixIndexT i = 0; i < 5; i++) {
3107  MatrixIndexT dimM = (Rand() % 10) + 10;
3108  MatrixIndexT dimO = (Rand() % 10) + 10;
3109 
3110  SpMatrix<Real> Q1(dimM), Q2(dimM), P1(dimO), P2(dimO);
3111  RandPosdefSpMatrix(dimM, &Q1);
3112  RandPosdefSpMatrix(dimM, &Q2);
3113  RandPosdefSpMatrix(dimO, &P1);
3114  RandPosdefSpMatrix(dimO, &P1);
3115 
3116  Matrix<Real> M(dimO, dimM), G(dimO, dimM);
3117  M.SetRandn();
3118  G.SetRandn();
3119  // InitRandNonsingular(&M);
3120  // InitRandNonsingular(&G);
3121 
3122  Matrix<Real> M2(M);
3123 
3124  SolverOptions opts;
3125  opts.optimize_delta = Rand() % 2;
3126  SolveDoubleQuadraticMatrixProblem(G, P1, P2, Q1, Q2, opts, &M2);
3127 
3128  {
3129  Real a1 = TraceMatMat(M2, G, kTrans),
3130  a2 = TraceMatSpMatSp(M2, kNoTrans, Q1, M2, kTrans, P1),
3131  a3 = TraceMatSpMatSp(M2, kNoTrans, Q2, M2, kTrans, P2),
3132  b1 = TraceMatMat(M, G, kTrans),
3133  b2 = TraceMatSpMatSp(M, kNoTrans, Q1, M, kTrans, P1),
3134  b3 = TraceMatSpMatSp(M, kNoTrans, Q2, M, kTrans, P2),
3135  a4 = a1 - 0.5 * a2 - 0.5 * a3,
3136  b4 = b1 - 0.5 * b2 - 0.5 * b3;
3137  KALDI_LOG << "a4 = " << a4 << ", b4 = " << b4;
3138  KALDI_ASSERT(a4 >= b4);
3139  } // Check objf not decreased.
3140  }
3141 }
3142 
3143 template<typename Real> static void UnitTestMax2() {
3144  for (MatrixIndexT i = 0; i < 2; i++) {
3145  MatrixIndexT M = 1 + Rand() % 10, N = 1 + Rand() % 10;
3146  Matrix<Real> A(M, N), B(M, N), C(M, N), D(M, N);
3147  A.SetRandn();
3148  B.SetRandn();
3149  for (MatrixIndexT r = 0; r < M; r++)
3150  for (MatrixIndexT c = 0; c < N; c++)
3151  C(r, c) = std::max(A(r, c), B(r, c));
3152  D.CopyFromMat(A);
3153  D.Max(B);
3154  AssertEqual(C, D);
3155  }
3156 }
3157 
3158 template<typename Real> static void UnitTestMaxAbsEig() {
3159  for (MatrixIndexT i = 0; i < 1; i++) {
3160  SpMatrix<Real> M(10);
3161  M.SetRandn();
3162  Matrix<Real> P(10, 10);
3163  Vector<Real> s(10);
3164  M.Eig(&s, (i == 0 ? static_cast<Matrix<Real>*>(NULL) : &P));
3165  Real max_eig = std::max(-s.Min(), s.Max());
3166  AssertEqual(max_eig, M.MaxAbsEig());
3167  }
3168 }
3169 
3170 template<typename Real> static void UnitTestLbfgs() {
3173  for (MatrixIndexT iter = 0; iter < 3; iter++) {
3174  bool minimize = (iter % 2 == 0);
3175  MatrixIndexT dim = 1 + Rand() % 30;
3176  SpMatrix<Real> S(dim);
3177  RandPosdefSpMatrix(dim, &S);
3178  Vector<Real> v(dim);
3179  v.SetRandn();
3180  // Function will be f = exp(0.1 * [ x' v -0.5 x' S x ])
3181  // This is to maximize; we negate it when minimizing.
3182 
3183  //Vector<Real> hessian(dim);
3184  //hessian.CopyDiagFromSp(S);
3185 
3186  SpMatrix<Real> Sinv(S);
3187  Sinv.Invert();
3188  Vector<Real> x_opt(dim);
3189  x_opt.AddSpVec(1.0, Sinv, v, 0.0); // S^{-1} v-- the optimum.
3190 
3191  Vector<Real> init_x(dim);
3192  init_x.SetRandn();
3193 
3194  LbfgsOptions opts;
3195  opts.minimize = minimize; // This objf has a maximum, not a minimum.
3196  OptimizeLbfgs<Real> opt_lbfgs(init_x, opts);
3197  MatrixIndexT num_iters = 0;
3198  Real c = 0.01;
3199  Real sign = (minimize ? -1.0 : 1.0); // function has a maximum not minimum..
3200  while (opt_lbfgs.RecentStepLength() > 1.0e-04) {
3201  KALDI_VLOG(2) << "Last step length is " << opt_lbfgs.RecentStepLength();
3202  const VectorBase<Real> &x = opt_lbfgs.GetProposedValue();
3203  Real logf = VecVec(x, v) - 0.5 * VecSpVec(x, S, x);
3204  Vector<Real> dlogf_dx(v); // derivative of log(f) w.r.t. x.
3205  dlogf_dx.AddSpVec(-1.0, S, x, 1.0);
3206  KALDI_VLOG(2) << "Gradient magnitude is " << dlogf_dx.Norm(2.0);
3207  Real f = Exp(c * logf);
3208  Vector<Real> df_dx(dlogf_dx);
3209  df_dx.Scale(f * c); // comes from derivative of the exponential function.
3210  f *= sign;
3211  df_dx.Scale(sign);
3212  opt_lbfgs.DoStep(f, df_dx);
3213  num_iters++;
3214  }
3215  Vector<Real> x (opt_lbfgs.GetValue());
3216  Vector<Real> diff(x);
3217  diff.AddVec(-1.0, x_opt);
3218  KALDI_VLOG(2) << "L-BFGS finished after " << num_iters << " function evaluations.";
3219  /*
3220  if (sizeof(Real) == 8) {
3221  KALDI_ASSERT(diff.Norm(2.0) < 0.5);
3222  } else {
3223  KALDI_ASSERT(diff.Norm(2.0) < 2.0);
3224  } */
3225  }
3226  g_kaldi_verbose_level = temp;
3227 }
3228 
3229 
3230 template<typename Real> static void UnitTestLinearCgd() {
3231  for (int i = 0; i < 20 ; i++) {
3232  MatrixIndexT M = 1 + Rand() % 20;
3233 
3234  SpMatrix<Real> A(M);
3235  RandPosdefSpMatrix(M, &A);
3236  Vector<Real> x(M), b(M), b2(M);
3237 
3238  LinearCgdOptions opts;
3239  if (Rand() % 2 == 0)
3240  opts.max_iters = 1 + Rand() % 10;
3241  if (Rand() % 2 == 0)
3242  opts.max_error = 1.0; // note: an absolute, not relative, error.
3243 
3244  x.SetRandn();
3245 
3246  b.AddSpVec(1.0, A, x, 0.0);
3247  Vector<Real> x_e(M); // x_e means x_estimated.
3248  x_e.SetRandn();
3249 
3250  int32 iters = LinearCgd(opts, A, b, &x_e);
3251 
3252  b2.AddSpVec(1.0, A, x_e, 0.0);
3253 
3254  Vector<Real> residual_error(b);
3255  residual_error.AddVec(-1.0, b2);
3256 
3257  BaseFloat error = residual_error.Norm(2.0);
3258 
3259  if (iters >= M) {
3260  // should have converged fully.
3261  Real max_abs = A.MaxAbsEig();
3262  KALDI_LOG << "error = " << error << ", b norm is " << b.Norm(2.0)
3263  << ", A max-abs-eig is " << max_abs;
3264  KALDI_ASSERT(error < 1.0e-04 * b.Norm(2.0) * max_abs);
3265  } else {
3266  BaseFloat wiggle_room = 1.1;
3267  if (opts.max_iters >= 0) {
3268  KALDI_ASSERT(iters <= opts.max_iters);
3269  if (iters < opts.max_iters) {
3270  KALDI_ASSERT(error <= wiggle_room * opts.max_error);
3271  }
3272  } else {
3273  KALDI_ASSERT(error <= wiggle_room * opts.max_error);
3274  }
3275  }
3276  }
3277 }
3278 
3279 
3280 template<typename Real> static void UnitTestMaxMin() {
3281 
3282  MatrixIndexT M = 1 + Rand() % 10, N = 1 + Rand() % 10;
3283  {
3284  Vector<Real> v(N);
3285  v.SetRandn();
3286  Real min = 1.0e+10, max = -1.0e+10;
3287  for (MatrixIndexT i = 0; i< N; i++) {
3288  min = std::min(min, v(i));
3289  max = std::max(max, v(i));
3290  }
3291  AssertEqual(min, v.Min());
3292  AssertEqual(max, v.Max());
3293  }
3294  {
3295  SpMatrix<Real> S(N);
3296  S.SetRandn();
3297  Real min = 1.0e+10, max = -1.0e+10;
3298  for (MatrixIndexT i = 0; i< N; i++) {
3299  for (MatrixIndexT j = 0; j <= i; j++) {
3300  min = std::min(min, S(i, j));
3301  max = std::max(max, S(i, j));
3302  }
3303  }
3304  AssertEqual(min, S.Min());
3305  AssertEqual(max, S.Max());
3306  }
3307  {
3308  Matrix<Real> mat(M, N);
3309  mat.SetRandn();
3310  Real min = 1.0e+10, max = -1.0e+10;
3311  for (MatrixIndexT i = 0; i< M; i++) {
3312  for (MatrixIndexT j = 0; j < N; j++) {
3313  min = std::min(min, mat(i, j));
3314  max = std::max(max, mat(i, j));
3315  }
3316  }
3317  AssertEqual(min, mat.Min());
3318  AssertEqual(max, mat.Max());
3319  }
3320 }
3321 
3322 template<typename Real>
3323 static bool approx_equal(Real a, Real b) {
3324  return ( std::abs(a-b) <= 1.0e-03 * (std::abs(a)+std::abs(b)));
3325 }
3326 
3327 template<typename Real> static void UnitTestTrace() {
3328 
3329  for (MatrixIndexT i = 0;i < 5;i++) {
3330  MatrixIndexT dimM = 20 + Rand()%10, dimN = 20 + Rand()%10, dimO = 20 + Rand()%10, dimP = dimM;
3331  Matrix<Real> A(dimM, dimN), B(dimN, dimO), C(dimO, dimP);
3332  A.SetRandn(); B.SetRandn(); C.SetRandn();
3333  Matrix<Real> AT(dimN, dimM), BT(dimO, dimN), CT(dimP, dimO);
3334  AT.CopyFromMat(A, kTrans); BT.CopyFromMat(B, kTrans); CT.CopyFromMat(C, kTrans);
3335 
3336  Matrix<Real> AB(dimM, dimO);
3337  AB.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0);
3338  Matrix<Real> BC(dimN, dimP);
3339  BC.AddMatMat(1.0, B, kNoTrans, C, kNoTrans, 0.0);
3340  Matrix<Real> ABC(dimM, dimP);
3341  ABC.AddMatMat(1.0, A, kNoTrans, BC, kNoTrans, 0.0);
3342 
3343  Real
3344  t1 = TraceMat(ABC),
3345  t2 = ABC.Trace(),
3346  t3 = TraceMatMat(A, BC),
3347  t4 = TraceMatMat(AT, BC, kTrans),
3348  t5 = TraceMatMat(BC, AT, kTrans),
3349  t6 = TraceMatMatMat(A, kNoTrans, B, kNoTrans, C, kNoTrans),
3350  t7 = TraceMatMatMat(AT, kTrans, B, kNoTrans, C, kNoTrans),
3351  t8 = TraceMatMatMat(AT, kTrans, BT, kTrans, C, kNoTrans),
3352  t9 = TraceMatMatMat(AT, kTrans, BT, kTrans, CT, kTrans);
3353 
3354  Matrix<Real> ABC1(dimM, dimP); // tests AddMatMatMat.
3355  ABC1.AddMatMatMat(1.0, A, kNoTrans, B, kNoTrans, C, kNoTrans, 0.0);
3356  AssertEqual(ABC, ABC1);
3357 
3358  Matrix<Real> ABC2(dimM, dimP); // tests AddMatMatMat.
3359  ABC2.AddMatMatMat(0.25, A, kNoTrans, B, kNoTrans, C, kNoTrans, 0.0);
3360  ABC2.AddMatMatMat(0.25, AT, kTrans, B, kNoTrans, C, kNoTrans, 2.0); // the extra 1.0 means another 0.25.
3361  ABC2.AddMatMatMat(0.125, A, kNoTrans, BT, kTrans, C, kNoTrans, 1.0);
3362  ABC2.AddMatMatMat(0.125, A, kNoTrans, B, kNoTrans, CT, kTrans, 1.0);
3363  AssertEqual(ABC, ABC2);
3364 
3365  Real tol = 0.001;
3366  KALDI_ASSERT((std::abs(t1-t2) < tol) && (std::abs(t2-t3) < tol) && (std::abs(t3-t4) < tol)
3367  && (std::abs(t4-t5) < tol) && (std::abs(t5-t6) < tol) && (std::abs(t6-t7) < tol)
3368  && (std::abs(t7-t8) < tol) && (std::abs(t8-t9) < tol));
3369  }
3370 
3371  for (MatrixIndexT i = 0;i < 5;i++) {
3372  MatrixIndexT dimM = 20 + Rand()%10, dimN = 20 + Rand()%10;
3373  SpMatrix<Real> S(dimM), T(dimN);
3374  S.SetRandn(); T.SetRandn();
3375  Matrix<Real> M(dimM, dimN), O(dimM, dimN);
3376  M.SetRandn(); O.SetRandn();
3377  Matrix<Real> sM(S), tM(T);
3378 
3379  Real x1 = TraceMatMat(tM, tM);
3380  Real x2 = TraceSpMat(T, tM);
3381  KALDI_ASSERT(approx_equal(x1, x2) || fabs(x1-x2) < 0.1);
3382 
3383  Real t1 = TraceMatMatMat(M, kNoTrans, tM, kNoTrans, M, kTrans);
3384  Real t2 = TraceMatSpMat(M, kNoTrans, T, M, kTrans);
3385  KALDI_ASSERT(approx_equal(t1, t2) || fabs(t1-12) < 0.1);
3386 
3387  Real u1 = TraceMatSpMatSp(M, kNoTrans, T, O, kTrans, S);
3388  Real u2 = TraceMatMatMatMat(M, kNoTrans, tM, kNoTrans, O, kTrans, sM, kNoTrans);
3389  KALDI_ASSERT(approx_equal(u1, u2) || fabs(u1-u2) < 0.1);
3390  }
3391 
3392 }
3393 
3394 
3395 template<typename Real> static void UnitTestComplexFt() {
3396 
3397  // Make sure it inverts properly.
3398  for (MatrixIndexT d = 0; d < 10; d++) {
3399  MatrixIndexT N = Rand() % 100, twoN = 2*N;
3400  Vector<Real> v(twoN), w(twoN), x(twoN);
3401  v.SetRandn();
3402  ComplexFt(v, &w, true);
3403  ComplexFt(w, &x, false);
3404  if (N>0) x.Scale(1.0/static_cast<Real>(N));
3405  AssertEqual(v, x);
3406  }
3407 }
3408 
3409 template<typename Real> static void UnitTestDct() {
3410 
3411  // Check that DCT matrix is orthogonal (i.e. M^T M = I);
3412  for (MatrixIndexT i = 0; i < 10; i++) {
3413  MatrixIndexT N = 1 + Rand() % 10;
3414  Matrix<Real> M(N, N);
3415  ComputeDctMatrix(&M);
3416  Matrix<Real> I(N, N);
3417  I.AddMatMat(1.0, M, kTrans, M, kNoTrans, 0.0);
3418  KALDI_ASSERT(I.IsUnit());
3419  }
3420 
3421 }
3422 template<typename Real> static void UnitTestComplexFft() {
3423 
3424  // Make sure it inverts properly.
3425  for (MatrixIndexT N_ = 0; N_ < 100; N_+=3) {
3426  MatrixIndexT N = N_;
3427  if (N>=95) {
3428  N = ( Rand() % 150);
3429  N = N*N; // big number.
3430  }
3431 
3432  MatrixIndexT twoN = 2*N;
3433  Vector<Real> v(twoN), w_base(twoN), w_alg(twoN), x_base(twoN), x_alg(twoN);
3434 
3435  v.SetRandn();
3436 
3437  if (N< 100) ComplexFt(v, &w_base, true);
3438  w_alg.CopyFromVec(v);
3439  ComplexFft(&w_alg, true);
3440  if (N< 100) AssertEqual(w_base, w_alg, 0.01*N);
3441 
3442  if (N< 100) ComplexFt(w_base, &x_base, false);
3443  x_alg.CopyFromVec(w_alg);
3444  ComplexFft(&x_alg, false);
3445 
3446  if (N< 100) AssertEqual(x_base, x_alg, 0.01*N);
3447  x_alg.Scale(1.0/N);
3448  AssertEqual(v, x_alg, 0.001*N);
3449  }
3450 }
3451 
3452 
3453 template<typename Real> static void UnitTestSplitRadixComplexFft() {
3454 
3455  // Make sure it inverts properly.
3456  for (MatrixIndexT N_ = 0; N_ < 30; N_+=3) {
3457  MatrixIndexT logn = 1 + Rand() % 10;
3458  MatrixIndexT N = 1 << logn;
3459 
3460  MatrixIndexT twoN = 2*N;
3461  std::vector<Real> temp_buffer;
3462  SplitRadixComplexFft<Real> srfft(N), srfft2(srfft);
3463  for (MatrixIndexT p = 0; p < 3; p++) {
3464  Vector<Real> v(twoN), w_base(twoN), w_alg(twoN), x_base(twoN), x_alg(twoN);
3465 
3466  v.SetRandn();
3467 
3468  if (N< 100) ComplexFt(v, &w_base, true);
3469  w_alg.CopyFromVec(v);
3470 
3471  if (Rand() % 2 == 0)
3472  srfft.Compute(w_alg.Data(), true);
3473  else
3474  srfft2.Compute(w_alg.Data(), true, &temp_buffer);
3475 
3476  if (N< 100) AssertEqual(w_base, w_alg, 0.01*N);
3477 
3478  if (N< 100) ComplexFt(w_base, &x_base, false);
3479  x_alg.CopyFromVec(w_alg);
3480  srfft.Compute(x_alg.Data(), false);
3481 
3482  if (N< 100) AssertEqual(x_base, x_alg, 0.01*N);
3483  x_alg.Scale(1.0/N);
3484  AssertEqual(v, x_alg, 0.001*N);
3485  }
3486  }
3487 }
3488 
3489 
3490 
3491 template<typename Real> static void UnitTestTranspose() {
3492 
3493  Matrix<Real> M(Rand() % 5 + 1, Rand() % 10 + 1);
3494  M.SetRandn();
3495  Matrix<Real> N(M, kTrans);
3496  N.Transpose();
3497  AssertEqual(M, N);
3498 }
3499 
3500 template<typename Real> static void UnitTestAddVecToRows() {
3501  std::vector<Real> sizes;
3502  sizes.push_back(16);
3503  sizes.push_back(128);
3504  for (int i = 0; i < 2; i++) {
3505  MatrixIndexT dimM = sizes[i] + Rand() % 10, dimN = sizes[i] + Rand() % 10;
3506  Matrix<Real> M(dimM, dimN);
3507  M.SetRandn();
3508  Vector<float> v(M.NumCols());
3509  v.SetRandn();
3510  Matrix<Real> N(M);
3511  Vector<float> ones(M.NumRows());
3512  ones.Set(1.0);
3513  M.AddVecToRows(0.5, v);
3514  N.AddVecVec(0.5, ones, v);
3515  AssertEqual(M, N);
3516  }
3517 }
3518 
3519 template<typename Real> static void UnitTestAddVec2Sp() {
3520  for (int32 i = 0; i < 10; i++) {
3521  int32 dim = Rand() % 5;
3522  SpMatrix<Real> S(dim);
3523  S.SetRandn();
3524  Vector<Real> v(dim);
3525  v.SetRandn();
3526  Matrix<Real> M(dim, dim);
3527  M.CopyDiagFromVec(v);
3528 
3529  SpMatrix<Real> T1(dim);
3530  T1.SetRandn();
3531  SpMatrix<Real> T2(T1);
3532  Real alpha = 0.33, beta = 4.5;
3533  T1.AddVec2Sp(alpha, v, S, beta);
3534  T2.AddMat2Sp(alpha, M, kNoTrans, S, beta);
3535  AssertEqual(T1, T2);
3536  }
3537 }
3538 
3539 
3540 template<typename Real> static void UnitTestAddVecToCols() {
3541  std::vector<Real> sizes;
3542  sizes.push_back(16);
3543  sizes.push_back(128);
3544  for (int i = 0; i < 2; i++) {
3545  MatrixIndexT dimM = sizes[i] + Rand() % 10, dimN = sizes[i] + Rand() % 10;
3546  Matrix<Real> M(dimM, dimN);
3547  M.SetRandn();
3548  Vector<float> v(M.NumRows());
3549  v.SetRandn();
3550  Matrix<Real> N(M);
3551  Vector<float> ones(M.NumCols());
3552  ones.Set(1.0);
3553  M.AddVecToCols(0.5, v);
3554  N.AddVecVec(0.5, v, ones);
3555  AssertEqual(M, N);
3556  }
3557 }
3558 
3559 template<typename Real> static void UnitTestComplexFft2() {
3560 
3561  // Make sure it inverts properly.
3562  for (MatrixIndexT pos = 0; pos < 10; pos++) {
3563  for (MatrixIndexT N_ = 2; N_ < 15; N_+=2) {
3564  if ( pos < N_) {
3565  MatrixIndexT N = N_;
3566  Vector<Real> v(N), vorig(N), v2(N);
3567  v(pos) = 1.0;
3568  vorig.CopyFromVec(v);
3569  // KALDI_LOG << "Original v:\n" << v;
3570  ComplexFft(&v, true);
3571  // KALDI_LOG << "one fft:\n" << v;
3572  ComplexFt(vorig, &v2, true);
3573  // KALDI_LOG << "one fft[baseline]:\n" << v2;
3574  if (!ApproxEqual(v, v2) ) {
3575  ComplexFft(&vorig, true);
3576  KALDI_ASSERT(0);
3577  }
3578  ComplexFft(&v, false);
3579  // KALDI_LOG << "one more:\n" << v;
3580  v.Scale(1.0/(N/2));
3581  if (!ApproxEqual(v, vorig)) {
3582  ComplexFft(&vorig, true);
3583  KALDI_ASSERT(0);
3584  }// AssertEqual(v, vorig);
3585  }
3586  }
3587  }
3588 }
3589 
3590 
3591 template<typename Real> static void UnitTestSplitRadixComplexFft2() {
3592 
3593  // Make sure it inverts properly.
3594  for (MatrixIndexT p = 0; p < 30; p++) {
3595  MatrixIndexT logn = 1 + Rand() % 10;
3596  MatrixIndexT N = 1 << logn;
3597  SplitRadixComplexFft<Real> srfft(N);
3598  for (MatrixIndexT q = 0; q < 3; q++) {
3599  Vector<Real> v(N*2), vorig(N*2);
3600  v.SetRandn();
3601  vorig.CopyFromVec(v);
3602  srfft.Compute(v.Data(), true); // forward
3603  srfft.Compute(v.Data(), false); // backward
3604  v.Scale(1.0/N);
3605  KALDI_ASSERT(ApproxEqual(v, vorig));
3606  }
3607  }
3608 }
3609 
3610 
3611 template<typename Real> static void UnitTestRealFft() {
3612 
3613  // First, test RealFftInefficient.
3614  for (MatrixIndexT N_ = 2; N_ < 100; N_ += 6) {
3615  MatrixIndexT N = N_;
3616  if (N >90) N *= Rand() % 60;
3617  Vector<Real> v(N), w(N), x(N), y(N);
3618  v.SetRandn();
3619  w.CopyFromVec(v);
3620  RealFftInefficient(&w, true);
3621  y.CopyFromVec(v);
3622  RealFft(&y, true); // test efficient one.
3623  // KALDI_LOG <<"v = "<<v;
3624  // KALDI_LOG << "Inefficient real fft of v is: "<< w;
3625  // KALDI_LOG << "Efficient real fft of v is: "<< y;
3626  AssertEqual(w, y, 0.01*N);
3627  x.CopyFromVec(w);
3628  RealFftInefficient(&x, false);
3629  RealFft(&y, false);
3630  // KALDI_LOG << "Inefficient real fft of v twice is: "<< x;
3631  if (N != 0) x.Scale(1.0/N);
3632  if (N != 0) y.Scale(1.0/N);
3633  AssertEqual(v, x, 0.001*N);
3634  AssertEqual(v, y, 0.001*N); // ?
3635  }
3636 }
3637 
3638 
3639 template<typename Real> static void UnitTestSplitRadixRealFft() {
3640 
3641  for (MatrixIndexT p = 0; p < 30; p++) {
3642  MatrixIndexT logn = 2 + Rand() % 9,
3643  N = 1 << logn;
3644 
3645  SplitRadixRealFft<Real> srfft(N), srfft2(srfft);
3646  std::vector<Real> temp_buffer;
3647  for (MatrixIndexT q = 0; q < 3; q++) {
3648  Vector<Real> v(N), w(N), x(N), y(N);
3649  v.SetRandn();
3650  w.CopyFromVec(v);
3651  RealFftInefficient(&w, true);
3652  y.CopyFromVec(v);
3653  if (Rand() % 2 == 0)
3654  srfft.Compute(y.Data(), true);
3655  else
3656  srfft2.Compute(y.Data(), true, &temp_buffer);
3657 
3658  // KALDI_LOG <<"v = "<<v;
3659  // KALDI_LOG << "Inefficient real fft of v is: "<< w;
3660  // KALDI_LOG << "Efficient real fft of v is: "<< y;
3661  AssertEqual(w, y, 0.01*N);
3662  x.CopyFromVec(w);
3663  RealFftInefficient(&x, false);
3664  srfft.Compute(y.Data(), false);
3665  // KALDI_LOG << "Inefficient real fft of v twice is: "<< x;
3666  x.Scale(1.0/N);
3667  y.Scale(1.0/N);
3668  AssertEqual(v, x, 0.001*N);
3669  AssertEqual(v, y, 0.001*N); // ?
3670  }
3671  }
3672 }
3673 
3674 
3675 
3676 template<typename Real> static void UnitTestRealFftSpeed() {
3677 
3678  // First, test RealFftInefficient.
3679  KALDI_LOG << "starting. ";
3680  MatrixIndexT sz = 512; // fairly typical size.
3681  for (MatrixIndexT i = 0; i < 3000; i++) {
3682  if (i % 1000 == 0) KALDI_LOG << "done 1000 [ == ten seconds of speech]";
3683  Vector<Real> v(sz);
3684  RealFft(&v, true);
3685  }
3686 }
3687 
3688 template<typename Real> static void UnitTestSplitRadixRealFftSpeed() {
3689  KALDI_LOG << "starting. ";
3690  MatrixIndexT sz = 512; // fairly typical size.
3691  SplitRadixRealFft<Real> srfft(sz);
3692  for (MatrixIndexT i = 0; i < 6000; i++) {
3693  if (i % 1000 == 0)
3694  KALDI_LOG << "done 1000 [ == ten seconds of speech, split-radix]";
3695  Vector<Real> v(sz);
3696  srfft.Compute(v.Data(), true);
3697  }
3698 }
3699 
3700 template<typename Real>
3702  // This tests a not-really-public function that's used in Matrix::Power().
3703 
3704  for (MatrixIndexT i = 0; i < 10; i++) {
3705  Real power = RandGauss();
3706  Real x = 2.0, y = 0.0;
3707  bool ans = AttemptComplexPower(&x, &y, power);
3708  KALDI_ASSERT(ans);
3709  AssertEqual(std::pow(static_cast<Real>(2.0), power), x);
3710  AssertEqual(y, 0.0);
3711  }
3712  {
3713  Real x, y;
3714  x = 0.5; y = -0.3;
3715  bool ans = AttemptComplexPower(&x, &y, static_cast<Real>(2.21));
3716  KALDI_ASSERT(ans);
3717  ans = AttemptComplexPower(&x, &y, static_cast<Real>(1.0/2.21));
3718  KALDI_ASSERT(ans);
3719  AssertEqual(x, 0.5);
3720  AssertEqual(y, -0.3);
3721  }
3722  {
3723  Real x, y;
3724  x = 0.5; y = -0.3;
3725  bool ans = AttemptComplexPower(&x, &y, static_cast<Real>(2.0));
3726  KALDI_ASSERT(ans);
3727  AssertEqual(x, 0.5*0.5 - 0.3*0.3);
3728  AssertEqual(y, -0.3*0.5*2.0);
3729  }
3730 
3731  {
3732  Real x, y;
3733  x = 1.0/std::sqrt(2.0); y = -1.0/std::sqrt(2.0);
3734  bool ans = AttemptComplexPower(&x, &y, static_cast<Real>(-1.0));
3735  KALDI_ASSERT(ans);
3736  AssertEqual(x, 1.0/std::sqrt(2.0));
3737  AssertEqual(y, 1.0/std::sqrt(2.0));
3738  }
3739 
3740  {
3741  Real x, y;
3742  x = 0.0; y = 0.0;
3743  bool ans = AttemptComplexPower(&x, &y, static_cast<Real>(-2.0));
3744  KALDI_ASSERT(!ans); // zero; negative pow.
3745  }
3746  {
3747  Real x, y;
3748  x = -2.0; y = 0.0;
3749  bool ans = AttemptComplexPower(&x, &y, static_cast<Real>(1.5));
3750  KALDI_ASSERT(!ans); // negative real case
3751  }
3752 }
3753 template<typename Real>
3755 
3756  for (MatrixIndexT iter = 0; iter < 30; iter++) {
3757  MatrixIndexT dimM = 1 + Rand() % 20;
3758  Matrix<Real> M(dimM, dimM);
3759  M.SetRandn();
3760 
3761  Matrix<Real> MM(dimM, dimM);
3762  MM.AddMatMat(1.0, M, kNoTrans, M, kNoTrans, 0.0); // MM = M M.
3763  Matrix<Real> MMMM(dimM, dimM);
3764  MMMM.AddMatMat(1.0, MM, kNoTrans, MM, kNoTrans, 0.0);
3765 
3766  Matrix<Real> MM2(MM);
3767  bool b = MM2.Power(1.0);
3768  KALDI_ASSERT(b);
3769  AssertEqual(MM2, MM);
3770  Matrix<Real> MMMM2(MM);
3771  b = MMMM2.Power(2.0);
3772  KALDI_ASSERT(b);
3773  AssertEqual(MMMM2, MMMM);
3774  }
3775  for (MatrixIndexT iter = 0; iter < 30; iter++) {
3776  MatrixIndexT dimM = 1 + Rand() % 20;
3777  Matrix<Real> M(dimM, dimM);
3778  InitRandNonsingular(&M);
3779 
3780  Matrix<Real> MM(dimM, dimM);
3781  MM.AddMatMat(1.0, M, kNoTrans, M, kNoTrans, 0.0); // MM = M M.
3782  // This ensures there are no real, negative eigenvalues.
3783 
3784  Matrix<Real> MMMM(dimM, dimM);
3785  MMMM.AddMatMat(1.0, MM, kNoTrans, MM, kNoTrans, 0.0);
3786 
3787  Matrix<Real> MM2(M);
3788  if (!MM2.Power(2.0)) { // possibly had negative eigenvalues
3789  KALDI_LOG << "Could not take matrix to power (not an error)";
3790  } else {
3791  AssertEqual(MM2, MM);
3792  }
3793  Matrix<Real> MMMM2(M);
3794  if (!MMMM2.Power(4.0)) { // possibly had negative eigenvalues
3795  KALDI_LOG << "Could not take matrix to power (not an error)";
3796  } else {
3797  AssertEqual(MMMM2, MMMM);
3798  }
3799  Matrix<Real> MMMM3(MM);
3800  if (!MMMM3.Power(2.0)) {
3801  KALDI_ERR << "Could not take matrix to power (should have been able to)";
3802  } else {
3803  AssertEqual(MMMM3, MMMM);
3804  }
3805 
3806  Matrix<Real> MM4(MM);
3807  if (!MM4.Power(-1.0))
3808  KALDI_ERR << "Could not take matrix to power (should have been able to)";
3809  MM4.Invert();
3810  AssertEqual(MM4, MM);
3811  }
3812 }
3813 
3815 
3816  Vector<float> v(5);
3817  v.SetRandn();
3818  Vector<double> w(5);
3819  w.SetRandn();
3820 
3821  Vector<float> wf(w);
3822 
3823  for (MatrixIndexT i = 0; i < 2; i++) {
3824  float f = 1.0;
3825  if (i == 0) f = 2.0;
3826 
3827  {
3828  Vector<float> sum1(5);
3829  Vector<double> sum2(5);
3830  Vector<float> sum3(5);
3831  sum1.AddVec(f, v); sum1.AddVec(f, w);
3832  sum2.AddVec(f, v); sum2.AddVec(f, w);
3833  sum3.AddVec(f, v); sum3.AddVec(f, wf);
3834  Vector<float> sum2b(sum2);
3835  AssertEqual(sum1, sum2b);
3836  AssertEqual(sum1, sum3);
3837  }
3838 
3839  {
3840  Vector<float> sum1(5);
3841  Vector<double> sum2(5);
3842  Vector<float> sum3(5);
3843  sum1.AddVec2(f, v); sum1.AddVec2(f, w);
3844  sum2.AddVec2(f, v); sum2.AddVec2(f, w);
3845  sum3.AddVec2(f, v); sum3.AddVec2(f, wf);
3846  Vector<float> sum2b(sum2);
3847  AssertEqual(sum1, sum2b);
3848  AssertEqual(sum1, sum3);
3849  }
3850  }
3851 }
3852 
3853 template<typename Real>
3854 static void UnitTestPca(bool full_test) {
3855  // We'll test that we can exactly reconstruct the vectors, if
3856  // the PCA dim is <= the "real" dim that the vectors live in.
3857  for (MatrixIndexT i = 0; i < 10; i++) {
3858  bool exact = i % 2 == 0;
3859  MatrixIndexT true_dim = (full_test ? 200 : 50) + Rand() % 5, //dim of subspace points live in
3860  feat_dim = true_dim + Rand() % 5, // dim of feature space
3861  num_points = true_dim + Rand() % 5, // number of training points.
3862  G = std::min(feat_dim,
3863  std::min(num_points,
3864  static_cast<MatrixIndexT>(true_dim + Rand() % 5)));
3865 
3866  Matrix<Real> Proj(feat_dim, true_dim);
3867  Proj.SetRandn();
3868  Matrix<Real> true_X(num_points, true_dim);
3869  true_X.SetRandn();
3870  Matrix<Real> X(num_points, feat_dim);
3871  X.AddMatMat(1.0, true_X, kNoTrans, Proj, kTrans, 0.0);
3872 
3873  Matrix<Real> U(G, feat_dim); // the basis
3874  Matrix<Real> A(num_points, G); // projection of points into the basis..
3875  ComputePca(X, &U, &A, true, exact);
3876  {
3877  SpMatrix<Real> I(G);
3878  I.AddMat2(1.0, U, kNoTrans, 0.0);
3879  KALDI_LOG << "Non-unit-ness of U is " << NonUnitness(I);
3880  KALDI_ASSERT(I.IsUnit(0.001));
3881  }
3882  Matrix<Real> X2(num_points, feat_dim);
3883  X2.AddMatMat(1.0, A, kNoTrans, U, kNoTrans, 0.0);
3884  // Check reproduction.
3885  KALDI_LOG << "A.Sum() " << A.Sum() << ", U.Sum() " << U.Sum();
3886  AssertEqual(X, X2, 0.01);
3887  // Check basis is orthogonal.
3888  Matrix<Real> tmp(G, G);
3889  tmp.AddMatMat(1.0, U, kNoTrans, U, kTrans, 0.0);
3890  KALDI_ASSERT(tmp.IsUnit(0.01));
3891  }
3892 }
3893 
3894 
3895 /* UnitTestPca2 test the same function, but it's more geared towards
3896  the 'inexact' method and for when you want less than the full number
3897  of PCA dimensions.
3898  */
3899 
3900 template<typename Real>
3901 static void UnitTestPca2(bool full_test) {
3902  for (MatrixIndexT i = 0; i < 5; i++) {
3903  // MatrixIndexT feat_dim = 600, num_points = 300;
3904  MatrixIndexT feat_dim = (full_test ? 600 : 100),
3905  num_points = (i%2 == 0 ? feat_dim * 2 : feat_dim / 2); // test
3906  // both branches of PCA code, inner + outer.
3907 
3908  Matrix<Real> X(num_points, feat_dim);
3909  X.SetRandn();
3910 
3911  MatrixIndexT pca_dim = 30;
3912 
3913  Matrix<Real> U(pca_dim, feat_dim); // rows PCA directions.
3914  bool print_eigs = true, exact = false;
3915  ComputePca(X, &U, static_cast<Matrix<Real>*>(NULL), print_eigs, exact);
3916 
3917  Real non_orth = NonOrthogonality(U, kNoTrans);
3918  KALDI_ASSERT(non_orth < 0.001);
3919  KALDI_LOG << "Non-orthogonality of U is " << non_orth;
3920  Matrix<Real> U2(pca_dim, feat_dim);
3921 
3922  {
3923  SpMatrix<Real> Scatter(feat_dim);
3924  Scatter.AddMat2(1.0, X, kTrans, 0.0);
3925  Matrix<Real> V(feat_dim, feat_dim);
3926  Vector<Real> l(feat_dim);
3927  Scatter.Eig(&l, &V); // cols of V are eigenvectors.
3928  SortSvd(&l, &V); // Get top dims.
3929  U2.CopyFromMat(SubMatrix<Real>(V, 0, feat_dim, 0, pca_dim), kTrans);
3930  Real non_orth = NonOrthogonality(U2, kNoTrans);
3931  KALDI_ASSERT(non_orth < 0.001);
3932  KALDI_LOG << "Non-orthogonality of U2 is " << non_orth;
3933 
3934  SpMatrix<Real> ScatterProjU(pca_dim), ScatterProjU2(pca_dim);
3935  ScatterProjU.AddMat2Sp(1.0, U, kNoTrans, Scatter, 0.0);
3936  ScatterProjU2.AddMat2Sp(1.0, U2, kNoTrans, Scatter, 0.0);
3937  KALDI_LOG << "Non-diagonality of proj with U is "
3938  << NonDiagonalness(ScatterProjU);
3939  KALDI_LOG << "Non-diagonality of proj with U2 is "
3940  << NonDiagonalness(ScatterProjU2);
3941  KALDI_ASSERT(ScatterProjU.IsDiagonal(0.01)); // Algorithm is statistical,
3942  // so it ends up being less accurate.
3943  KALDI_ASSERT(ScatterProjU2.IsDiagonal());
3944  KALDI_LOG << "Trace proj with U is " << ScatterProjU.Trace()
3945  << " with U2 is " << ScatterProjU2.Trace();
3946  AssertEqual(ScatterProjU.Trace(), ScatterProjU2.Trace(), 0.1);// Algorithm is
3947  // statistical so give it some leeway.
3948  }
3949  }
3950 }
3951 
3952 template<typename Real>
3953 static void UnitTestSvdSpeed() {
3954  std::vector<MatrixIndexT> sizes;
3955  sizes.push_back(100);
3956  sizes.push_back(150);
3957  sizes.push_back(200);
3958  sizes.push_back(300);
3959  sizes.push_back(500);
3960  sizes.push_back(750);
3961  sizes.push_back(1000);
3962  sizes.push_back(2000);
3963  time_t start, end;
3964  for (size_t i = 0; i < sizes.size(); i++) {
3965  MatrixIndexT size = sizes[i];
3966  {
3967  start = time(NULL);
3968  SpMatrix<Real> S(size);
3969  Vector<Real> l(size);
3970  S.Eig(&l);
3971  end = time(NULL);
3972  double diff = difftime(end, start);
3973  KALDI_LOG << "For size " << size << ", Eig without eigenvectors took " << diff
3974  << " seconds.";
3975  }
3976  {
3977  start = time(NULL);
3978  SpMatrix<Real> S(size);
3979  S.SetRandn();
3980  Vector<Real> l(size);
3981  Matrix<Real> P(size, size);
3982  S.Eig(&l, &P);
3983  end = time(NULL);
3984  double diff = difftime(end, start);
3985  KALDI_LOG << "For size " << size << ", Eig with eigenvectors took " << diff
3986  << " seconds.";
3987  }
3988  {
3989  start = time(NULL);
3990  Matrix<Real> M(size, size);
3991  M.SetRandn();
3992  Vector<Real> l(size);
3993  M.Svd(&l, NULL, NULL);
3994  end = time(NULL);
3995  double diff = difftime(end, start);
3996  KALDI_LOG << "For size " << size << ", SVD without eigenvectors took " << diff
3997  << " seconds.";
3998  }
3999  {
4000  start = time(NULL);
4001  Matrix<Real> M(size, size), U(size, size), V(size, size);
4002  M.SetRandn();
4003  Vector<Real> l(size);
4004  M.Svd(&l, &U, &V);
4005  end = time(NULL);
4006  double diff = difftime(end, start);
4007  KALDI_LOG << "For size " << size << ", SVD with eigenvectors took " << diff
4008  << " seconds.";
4009  }
4010  }
4011 }
4012 
4013 template<typename Real> static void UnitTestCompressedMatrix2() {
4014  // These are some new tests added after we add the capability to
4015  // specify the compression type.
4016 
4017  for (int32 i = 0; i < 10; i++) {
4018  // test that the kTwoByteSignedInteger method works.
4019  int32 num_rows = RandInt(1, 5), num_cols = RandInt(1, 10);
4020  Matrix<Real> mat(num_rows, num_cols);
4021  for (int32 j = 0; j < num_rows; j++) {
4022  for (int32 k = 0; k < num_cols; k++) {
4023  mat(j, k) = RandInt(-32768, 32767);
4024  }
4025  }
4027 
4028  Matrix<Real> mat2(cmat);
4029 
4030  // Check that they are exactly equal. These integers should all be
4031  // exactly representable, and exactly reconstructed after compression.
4032  KALDI_ASSERT(ApproxEqual(mat, mat2, Real(0.0)));
4033  }
4034 
4035 
4036  for (int32 i = 0; i < 10; i++) {
4037  // test that the kOneByteZeroOne compression method works.
4038  int32 num_rows = RandInt(1, 5), num_cols = RandInt(1, 10);
4039  Matrix<Real> mat(num_rows, num_cols);
4040  for (int32 j = 0; j < num_rows; j++) {
4041  for (int32 k = 0; k < num_cols; k++) {
4042  mat(j, k) = RandInt(0, 255) / 255.0;
4043  }
4044  }
4045  CompressedMatrix cmat(mat, kOneByteZeroOne);
4046 
4047  Matrix<Real> mat2(cmat);
4048 
4049  // Check that they are almost exactly equal. (It's not 100% exact because
4050  // 1.0 / 255.0 is not exactly representable.
4051  KALDI_ASSERT(ApproxEqual(mat, mat2, Real(1.00001)));
4052  }
4053 
4054 }
4055 
4056 
4057 template<typename Real> static void UnitTestCompressedMatrix() {
4058  // This is the basic test.
4059 
4060  CompressedMatrix empty_cmat; // some tests on empty matrix
4061  KALDI_ASSERT(empty_cmat.NumRows() == 0);
4062  KALDI_ASSERT(empty_cmat.NumCols() == 0);
4063 
4064  // could set num_tot to 10000 for more thorough testing.
4065  MatrixIndexT num_failure = 0, num_tot = 1000, max_failure = 1;
4066  for (MatrixIndexT n = 0; n < num_tot; n++) {
4067  MatrixIndexT num_rows = Rand() % 20, num_cols = Rand() % 15;
4068  if (num_rows * num_cols == 0) {
4069  num_rows = 0;
4070  num_cols = 0;
4071  }
4072  if (rand() % 2 == 0 && num_cols != 0) {
4073  // smaller matrices are more likely to have problems.
4074  num_cols = 1 + Rand() % 3;
4075  }
4076  Matrix<Real> M(num_rows, num_cols);
4077  if (Rand() % 3 != 0) M.SetRandn();
4078  else {
4079  M.Add(RandGauss());
4080  }
4081  if (Rand() % 2 == 0 && num_rows != 0) { // set one row to all the same value,
4082  // which is one possible pathology.
4083  // Give it large dynamic range to increase chance that it
4084  // is the largest or smallest value in the matrix.
4085  M.Row(Rand() % num_rows).Set(RandGauss() * 4.0);
4086  }
4087  double rand_val = RandGauss() * 4.0;
4088  // set a bunch of elements to all one value: increases
4089  // chance of pathologies.
4090  MatrixIndexT modulus = 1 + Rand() % 5;
4091  for (MatrixIndexT r = 0; r < num_rows; r++)
4092  for (MatrixIndexT c = 0; c < num_cols; c++)
4093  if (Rand() % modulus != 0) M(r, c) = rand_val;
4094 
4095 
4096  CompressionMethod method;
4097  switch(RandInt(0, 3)) {
4098  case 0: method = kAutomaticMethod; break;
4099  case 1: method = kSpeechFeature; break;
4100  case 2: method = kTwoByteAuto; break;
4101  default: method = kOneByteAuto; break;
4102  }
4103 
4104  CompressedMatrix cmat(M, method);
4105  KALDI_ASSERT(cmat.NumRows() == num_rows);
4106  KALDI_ASSERT(cmat.NumCols() == num_cols);
4107 
4108  Matrix<Real> M2(cmat.NumRows(), cmat.NumCols());
4109  cmat.CopyToMat(&M2);
4110 
4111  Matrix<Real> diff(M2);
4112  diff.AddMat(-1.0, M);
4113 
4114  { // Check that when compressing a matrix that has already been compressed,
4115  // and uncompressing, we get the same answer if using the same compression
4116  // method.
4117  // ok, actually, we can't guarantee this, so just limit the number of
4118  // failures.
4119  CompressedMatrix cmat2(M2, method);
4120  Matrix<Real> M3(cmat.NumRows(), cmat.NumCols());
4121  cmat2.CopyToMat(&M3);
4122  if (!M2.ApproxEqual(M3, 1.0e-04)) {
4123  KALDI_LOG << "cmat is: ";
4124  cmat.Write(std::cout, false);
4125  KALDI_LOG << "cmat2 is: ";
4126  cmat2.Write(std::cout, false);
4127  KALDI_WARN << "Matrices differ " << M2 << " vs. " << M3 << ", M2 range is "
4128  << M2.Min() << " to " << M2.Max() << ", M3 range is "
4129  << M3.Min() << " to " << M3.Max();
4130  num_failure++;
4131  }
4132  }
4133 
4134  if (num_rows > 0) { // Check that the constructor accepting row and column offsets
4135  // etc., works.
4136  if (RandInt(0, 1) == 0) {
4137  // test the ability of the self-constructor to do row-padding. (used in
4138  // getting nnet3 examples without un-compressing and re-compressing the
4139  // data_.
4140  bool allow_row_padding = true;
4141  int32 row_offset = RandInt(-4, num_rows - 1),
4142  col_offset = RandInt(0, num_cols - 1),
4143  num_rows_sub = RandInt(1, 4 + num_rows - row_offset),
4144  num_cols_sub = RandInt(1, num_cols - col_offset);
4145  CompressedMatrix cmat_sub(cmat, row_offset, num_rows_sub,
4146  col_offset, num_cols_sub, allow_row_padding);
4147  Matrix<Real> M2_sub(num_rows_sub, num_cols_sub);
4148  for (int32 row = 0; row < num_rows_sub; row++) {
4149  int32 old_row = row + row_offset;
4150  if (old_row < 0) old_row = 0;
4151  else if (old_row >= num_rows) { old_row = num_rows - 1; }
4152  SubVector<Real> M2_sub_row(M2_sub, row),
4153  M2_row(M2, old_row),
4154  M2_row_part(M2_row, col_offset, num_cols_sub);
4155  M2_sub_row.CopyFromVec(M2_row_part);
4156  }
4157  Matrix<Real> M3_sub(cmat_sub);
4158  M3_sub.AddMat(-1.0, M2_sub);
4159  KALDI_ASSERT(M3_sub.FrobeniusNorm() / (num_rows_sub * num_cols_sub) <
4160  1.0e-03);
4161  } else {
4162  int32 row_offset = RandInt(0, num_rows - 1),
4163  col_offset = RandInt(0, num_cols - 1),
4164  num_rows_sub = RandInt(1, num_rows - row_offset),
4165  num_cols_sub = RandInt(1, num_cols - col_offset);
4166  CompressedMatrix cmat_sub(cmat, row_offset, num_rows_sub,
4167  col_offset, num_cols_sub);
4168  SubMatrix<Real> M2_sub(M2, row_offset, num_rows_sub,
4169  col_offset, num_cols_sub);
4170  Matrix<Real> M3_sub(cmat_sub);
4171  M3_sub.AddMat(-1.0, M2_sub);
4172  KALDI_ASSERT(M3_sub.FrobeniusNorm() / (num_rows_sub * num_cols_sub) <
4173  1.0e-03);
4174  }
4175  } else {
4176  CompressedMatrix cmat_sub(cmat, 0, 0, 0, 0);
4177  }
4178 
4179  // test CopyRowToVec
4180  for (MatrixIndexT i = 0; i < num_rows; i++) {
4181  Vector<Real> V(num_cols);
4182  cmat.CopyRowToVec(i, &V); // get row.
4183  for (MatrixIndexT k = 0; k < num_cols; k++) {
4184  AssertEqual(M2(i, k), V(k));
4185  }
4186  }
4187 
4188  // test CopyColToVec
4189  for (MatrixIndexT i = 0; i < num_cols; i++) {
4190  Vector<Real> V(num_rows);
4191  cmat.CopyColToVec(i, &V); // get col.
4192  for (MatrixIndexT k = 0;k < num_rows;k++) {
4193  AssertEqual(M2(k, i), V(k));
4194  }
4195  }
4196 
4197  //test of getting a submatrix
4198  if(num_rows != 0 && num_cols != 0){
4199  MatrixIndexT sub_row_offset = Rand() % num_rows,
4200  sub_col_offset = Rand() % num_cols;
4201  // to make sure we don't mod by zero
4202  MatrixIndexT num_subrows = Rand() % (num_rows-sub_row_offset),
4203  num_subcols = Rand() % (num_cols-sub_col_offset);
4204  if(num_subrows == 0 || num_subcols == 0){ // in case we randomized to
4205  // empty matrix, at least make it correct
4206  num_subrows = 0;
4207  num_subcols = 0;
4208  }
4209  Matrix<Real> Msub(num_subrows, num_subcols);
4210  cmat.CopyToMat(sub_row_offset, sub_col_offset, &Msub);
4211  for (MatrixIndexT i = 0; i < num_subrows; i++) {
4212  for (MatrixIndexT k = 0;k < num_subcols;k++) {
4213  AssertEqual(M2(i+sub_row_offset, k+sub_col_offset), Msub(i, k));
4214  }
4215  }
4216  }
4217 
4218  { // Check Scale() method for compressedMatrix.
4219  for (int32 t = 0; t < 10; t++) {
4220  float alpha = 0.1;
4221  MatrixIndexT num_rows = 4 + Rand() % 20,
4222  num_cols = 10 + Rand() % 50;
4223  Matrix<Real> M(num_rows, num_cols);
4224  M.SetRandn();
4225  CompressedMatrix cmat(M);
4226  Matrix<Real> scaled_comp_mat(num_rows, num_cols),
4227  scaled_mat(M);
4228  scaled_mat.Scale(alpha);
4229  cmat.Scale(alpha);
4230  cmat.CopyToMat(&scaled_comp_mat);
4231  AssertEqual(scaled_comp_mat, scaled_mat);
4232  }
4233  }
4234  if (n < 5) { // test I/O.
4235  bool binary = (n % 2 == 1);
4236  {
4237  std::ofstream outs("tmpf", std::ios_base::out |std::ios_base::binary);
4238  InitKaldiOutputStream(outs, binary);
4239  cmat.Write(outs, binary);
4240  }
4241  CompressedMatrix cmat2;
4242  {
4243  bool binary_in;
4244  std::ifstream ins("tmpf", std::ios_base::in | std::ios_base::binary);
4245  InitKaldiInputStream(ins, &binary_in);
4246  cmat2.Read(ins, binary_in);
4247  }
4248 #if 1
4249  { // check that compressed-matrix can be read as matrix.
4250  bool binary_in;
4251  std::ifstream ins("tmpf", std::ios_base::in | std::ios_base::binary);
4252  InitKaldiInputStream(ins, &binary_in);
4253  Matrix<Real> mat1;
4254  mat1.Read(ins, binary_in);
4255  Matrix<Real> mat2(cmat2);
4256  AssertEqual(mat1, mat2);
4257  }
4258 #endif
4259 
4260 
4261  { // check that matrix can be read as compressed-matrix.
4262  Matrix<Real> mat1(cmat);
4263  {
4264  std::ofstream outs("tmpf", std::ios_base::out |std::ios_base::binary);
4265  InitKaldiOutputStream(outs, binary);
4266  mat1.Write(outs, binary);
4267  }
4268  bool binary_in;
4269  std::ifstream ins("tmpf", std::ios_base::in | std::ios_base::binary);
4270  InitKaldiInputStream(ins, &binary_in);
4271  CompressedMatrix cmat2;
4272  cmat2.Read(ins, binary_in);
4273  Matrix<Real> mat2(cmat2);
4274  AssertEqual(mat1, mat2);
4275  }
4276 
4277 
4278  Matrix<Real> M3(cmat2.NumRows(), cmat2.NumCols());
4279  cmat2.CopyToMat(&M3);
4280  AssertEqual(M2, M3); // tests I/O of CompressedMatrix.
4281 
4282  CompressedMatrix cmat3(cmat2); // testing self-constructor, which
4283  // tests assignment operator.
4284  Matrix<Real> M4(cmat3.NumRows(), cmat3.NumCols());
4285  cmat3.CopyToMat(&M4);
4286  AssertEqual(M2, M4);
4287  }
4288  KALDI_LOG << "M = " << M;
4289  KALDI_LOG << "M2 = " << M2;
4290  double tot = M.FrobeniusNorm(), err = diff.FrobeniusNorm();
4291  KALDI_LOG << "Compressed matrix, tot = " << tot << ", diff = "
4292  << err;
4293  if (err > 0.015 * tot) {
4294  KALDI_WARN << "Failure in compressed-matrix test.";
4295  num_failure++;
4296  }
4297  }
4298  if (num_failure > max_failure)
4299  KALDI_ERR << "Too many failures in compressed matrix test " << num_failure
4300  << " > " << max_failure;
4301 
4302  unlink("tmpf");
4303 }
4304 
4305 template<typename Real> static void UnitTestGeneralMatrix() {
4306  // This is the basic test.
4307 
4308  GeneralMatrix empty_pmat; // some tests on empty matrix
4309  KALDI_ASSERT(empty_pmat.NumRows() == 0);
4310  KALDI_ASSERT(empty_pmat.NumCols() == 0);
4311 
4312  // could set num_tot to 10000 for more thorough testing.
4313  MatrixIndexT num_failure = 0, num_tot = 1000, max_failure = 1;
4314  for (MatrixIndexT n = 0; n < num_tot; n++) {
4315  MatrixIndexT num_rows = Rand() % 20, num_cols = Rand() % 15;
4316  if (num_rows * num_cols == 0) {
4317  num_rows = 0;
4318  num_cols = 0;
4319  }
4320  if (rand() % 2 == 0 && num_cols != 0) {
4321  // smaller matrices are more likely to have problems.
4322  num_cols = 1 + Rand() % 3;
4323  }
4324  Matrix<Real> M(num_rows, num_cols);
4325  if (Rand() % 3 != 0)
4326  M.SetRandn();
4327  else {
4328  M.Add(RandGauss());
4329  }
4330  if (Rand() % 2 == 0 && num_rows != 0) { // set one row to all the same value,
4331  // which is one possible pathology.
4332  // Give it large dynamic range to increase chance that it
4333  // is the largest or smallest value in the matrix.
4334  M.Row(Rand() % num_rows).Set(RandGauss() * 4.0);
4335  }
4336  double rand_val = RandGauss() * 4.0;
4337  // set a bunch of elements to all one value: increases
4338  // chance of pathologies.
4339  MatrixIndexT modulus = 1 + Rand() % 5;
4340  for (MatrixIndexT r = 0; r < num_rows; r++)
4341  for (MatrixIndexT c = 0; c < num_cols; c++)
4342  if (Rand() % modulus != 0) M(r, c) = rand_val;
4343 
4344  GeneralMatrix pmat(M);
4345  if (RandInt(0, 1) == 0)
4346  pmat.Compress();
4347 
4348  if (RandInt(0, 1) == 0) {
4349  SparseMatrix<BaseFloat> smat(num_rows, num_cols);
4350  smat.SetRandn(0.1);
4351  pmat.Clear();
4352  smat.CopyToMat(&M, kNoTrans);
4353  pmat.SwapSparseMatrix(&smat);
4354  }
4355 
4356  KALDI_ASSERT(pmat.NumRows() == num_rows);
4357  KALDI_ASSERT(pmat.NumCols() == num_cols);
4358  GeneralMatrix pmat2(pmat);
4359 
4360  Matrix<Real> M2(pmat2.NumRows(), pmat2.NumCols());
4361  pmat2.GetMatrix(&M2);
4362 
4363  Matrix<Real> diff(M2);
4364  diff.AddMat(-1.0, M);
4365 
4366  if (pmat2.NumRows() > 0) { // Test ExtractRowRangeWithPadding.
4367  int32 row_offset = RandInt(-3, pmat.NumRows() - 1),
4368  num_rows = RandInt(1, 3 + pmat.NumRows() - row_offset);
4369  GeneralMatrix pmat3;
4370  ExtractRowRangeWithPadding(pmat2, row_offset, num_rows, &pmat3);
4371 
4372  Matrix<Real> mat_A(num_rows, pmat.NumCols()),
4373  mat_B(num_rows, pmat.NumCols());
4374  pmat3.GetMatrix(&mat_A);
4375  for (int32 row_out = 0; row_out < num_rows; row_out++) {
4376  int32 row_in = row_out + row_offset;
4377  if (row_in < 0) row_in = 0;
4378  if (row_in >= pmat2.NumRows())
4379  row_in = pmat2.NumRows() - 1;
4380  SubVector<Real> vec_out(mat_B, row_out),
4381  vec_in(M2, row_in);
4382  vec_out.CopyFromVec(vec_in);
4383  }
4384  if (mat_A.NumRows() >= 8) {
4385  // there should be exact equality.
4386  AssertEqual(mat_A, mat_B, 0.0);
4387  } else {
4388  // If it was compressed matrix and num-rows selected is < 8 and the
4389  // format was the one with per-column headers, then we would have changed
4390  // the format to save memory, which is not completely lossless.
4391  mat_A.AddMat(-1.0, mat_B);
4392  KALDI_ASSERT(mat_A.FrobeniusNorm() /
4393  (mat_A.NumRows() * mat_A.NumCols()) < 0.001);
4394  }
4395  }
4396 
4397  if (n < 5) { // test I/O.
4398  bool binary = (n % 2 == 1);
4399  {
4400  std::ofstream outs("tmpf", std::ios_base::out |std::ios_base::binary);
4401  InitKaldiOutputStream(outs, binary);
4402  pmat.Write(outs, binary);
4403  }
4404  GeneralMatrix pmat3;
4405  {
4406  bool binary_in;
4407  std::ifstream ins("tmpf", std::ios_base::in | std::ios_base::binary);
4408  InitKaldiInputStream(ins, &binary_in);
4409  pmat3.Read(ins, binary_in);
4410  }
4411 
4412  Matrix<Real> M3(pmat3.NumRows(), pmat3.NumCols());
4413  pmat3.GetMatrix(&M3);
4414  AssertEqual(M, M3);
4415  }
4416 
4417  KALDI_LOG << "M = " << M;
4418  KALDI_LOG << "M2 = " << M2;
4419  double tot = M.FrobeniusNorm(), err = diff.FrobeniusNorm();
4420  KALDI_LOG << "Compressed matrix, tot = " << tot << ", diff = "
4421  << err;
4422  if (err > 0.015 * tot) {
4423  KALDI_WARN << "Failure in possibly compressed-matrix test.";
4424  num_failure++;
4425  }
4426  }
4427  if (num_failure > max_failure)
4428  KALDI_ERR << "Too many failures in possibly compressed matrix test " << num_failure
4429  << " > " << max_failure;
4430 
4431  unlink("tmpf");
4432 }
4433 
4434 template<typename Real>
4436  for (int32 i = 0; i < 30; i++) {
4437  MatrixIndexT num_rows = Rand() % 20, num_cols = Rand() % 30;
4438  if (num_rows * num_cols == 0) {
4439  // this test wouldn't work for empty matrices.
4440  num_rows++;
4441  num_cols++;
4442  }
4443  Matrix<Real> mat(num_rows, num_cols);
4444  mat.SetRandn();
4445  CompressedMatrix cmat(mat);
4446 
4447  MatrixIndexT row_offset = Rand() % num_rows, col_offset = Rand() % num_cols;
4448  MatrixIndexT sub_num_rows = Rand() % (num_rows - row_offset) + 1,
4449  sub_num_cols = Rand() % (num_cols - col_offset) + 1;
4450  KALDI_VLOG(3) << "Whole matrix size: " << num_rows << "," << num_cols;
4451  KALDI_VLOG(3) << "Sub-matrix size: " << sub_num_rows << "," << sub_num_cols
4452  << " with offsets " << row_offset << "," << col_offset;
4453  CompressedMatrix cmat2(cmat, row_offset, sub_num_rows, //take a subset of
4454  col_offset, sub_num_cols); // the compressed matrix
4455  Matrix<Real> mat2(sub_num_rows, sub_num_cols);
4456  cmat2.CopyToMat(&mat2); // uncompress the subset of the compressed matrix
4457 
4458  Matrix<Real> mat3(cmat); // uncompress the whole compressed matrix
4459  SubMatrix<Real> sub_mat(mat3, row_offset, sub_num_rows, col_offset, sub_num_cols);
4460  if(!sub_mat.ApproxEqual(mat2)) {
4461  KALDI_ERR << "Matrices differ " << sub_mat << " vs. " << mat2;
4462  }
4463  }
4464 }
4465 
4466 
4467 template<typename Real>
4468 static void UnitTestTridiag() {
4469  SpMatrix<Real> A(3);
4470  A(1,1) = 1.0;
4471  A(1, 2) = 1.0;
4473  A(0, 2) = 1.0;
4475  A(0, 2) = 0.0;
4476  A(0, 1) = 1.0;
4478 }
4479 
4480 template<typename Real>
4482  int32 N = 1 + Rand() % 10;
4483  Vector<Real> vec(N);
4484  for (int32 n = 0; n < N; n++)
4485  vec(n) = Rand() % 3;
4486  if (vec.Sum() == 0)
4487  vec(0) = 2.0;
4488  Real sum = vec.Sum();
4489  int32 num_samples = 100000;
4490  std::vector<int32> counts(N, 0);
4491  for (int32 i = 0; i < num_samples; i++)
4492  counts[vec.RandCategorical()]++;
4493  for (int32 n = 0; n < N; n++) {
4494  Real a = counts[n] / (1.0 * num_samples),
4495  b = vec(n) / sum;
4496  KALDI_LOG << "a = " << a << ", b = " << b;
4497  KALDI_ASSERT(fabs(a - b) <= 0.1); // pretty arbitrary. Will increase #samp if fails.
4498  }
4499 }
4500 
4501 
4502 template <class Real>
4503 static void UnitTestAddMatMatNans() {
4504  for (int32 i = 0; i < 200; i++) {
4505  int32 num_rows = RandInt(1, 256), mid = RandInt(1, 256), num_cols = RandInt(1, 256);
4506  Matrix<Real> mat1(num_rows, mid), mat2(mid, num_cols), prod(num_rows, num_cols);
4507  PlaceNansInGaps(&mat1);
4508  PlaceNansInGaps(&mat2);
4509  prod.AddMatMat(1.0, mat1, kNoTrans, mat2, kNoTrans, 0.0);
4510  // make sure the nan's don't propagate.
4511  KALDI_ASSERT(prod.Sum() == 0.0 &&
4512  "The BLAS library that you are linking against has an issue that might "
4513  "cause problems later on.");
4514  }
4515 }
4516 
4517 template<class Real>
4518 static void UnitTestTopEigs() {
4519  for (MatrixIndexT i = 0; i < 2; i++) {
4520  // Previously tested with this but takes too long.
4521  MatrixIndexT dim = 400, num_eigs = 100;
4522  SpMatrix<Real> mat(dim);
4523  for (MatrixIndexT i = 0; i < dim; i++)
4524  for (MatrixIndexT j = 0; j <= i; j++)
4525  mat(i, j) = RandGauss();
4526 
4527  Matrix<Real> P(dim, num_eigs);
4528  Vector<Real> s(num_eigs);
4529  mat.TopEigs(&s, &P);
4530  { // P should have orthogonal columns. Check this.
4531  SpMatrix<Real> S(num_eigs);
4532  S.AddMat2(1.0, P, kTrans, 0.0);
4533  KALDI_LOG << "Non-unit-ness of S is " << NonUnitness(S);
4534  KALDI_ASSERT(S.IsUnit(1.0e-04));
4535  }
4536  // Note: we call the matrix "mat" by the name "S" below.
4537  Matrix<Real> SP(dim, num_eigs); // diag of P^T SP should be eigs.
4538  SP.AddSpMat(1.0, mat, P, kNoTrans, 0.0);
4539  Matrix<Real> PSP(num_eigs, num_eigs);
4540  PSP.AddMatMat(1.0, P, kTrans, SP, kNoTrans, 0.0);
4541  Vector<Real> s2(num_eigs);
4542  s2.CopyDiagFromMat(PSP);
4543  AssertEqual(s, s2);
4544  // Now check that eigs are close to real top eigs.
4545  {
4546  Matrix<Real> fullP(dim, dim);
4547  Vector<Real> fulls(dim);
4548  mat.Eig(&fulls, &fullP);
4549  KALDI_LOG << "Approximate eigs are " << s;
4550  // find sum of largest-abs-value eigs.
4551  fulls.ApplyAbs();
4552  std::sort(fulls.Data(), fulls.Data() + dim);
4553  SubVector<Real> tmp(fulls, dim - num_eigs, num_eigs);
4554  KALDI_LOG << "abs(real eigs) are " << tmp;
4555  BaseFloat real_sum = tmp.Sum();
4556  s.ApplyAbs();
4557  BaseFloat approx_sum = s.Sum();
4558  KALDI_LOG << "real sum is " << real_sum << ", approx_sum = " << approx_sum;
4559  }
4560  }
4561 }
4562 
4563 template<typename Real> static void UnitTestTriVecSolver() {
4564  for (MatrixIndexT iter = 0; iter < 100; iter++) {
4565  int32 dim = 1 + Rand() % 20;
4566  Vector<Real> b(dim);
4567  b.SetRandn();
4568  TpMatrix<Real> T(dim);
4569  T.SetRandn();
4570 
4571  bool bad = false;
4572  for (int32 i = 0; i < dim; i++) {
4573  if (fabs(T(i, i)) < 0.2)
4574  bad = true;
4575  }
4576  if (bad) {
4577  // Test may fail due to almost-singular matrix.
4578  continue;
4579  }
4580 
4581  Vector<Real> x(b);
4582  MatrixTransposeType trans = (iter % 2 == 0 ? kTrans : kNoTrans);
4583  x.Solve(T, trans); // solve for T x = b
4584  Vector<Real> b2(dim);
4585  b2.AddTpVec((Real)1.0, T, trans, x, (Real)0.0);
4586  KALDI_LOG << "b is " << b << ", b2 is " << b2;
4587  AssertEqual(b, b2, 0.01);
4588  }
4589 }
4590 
4591 
4592 template<typename Real> static void MatrixUnitTest(bool full_test) {
4593  UnitTestLinearCgd<Real>();
4594  UnitTestGeneralMatrix<BaseFloat>();
4595  UnitTestTridiagonalize<Real>();
4596  UnitTestTridiagonalizeAndQr<Real>();
4597  UnitTestAddMatSmat<Real>();
4598  UnitTestFloorChol<Real>();
4599  UnitTestFloorUnit<Real>();
4600  UnitTestAddMat2Sp<Real>();
4601  UnitTestLbfgs<Real>();
4602  // UnitTestSvdBad<Real>(); // test bug in Jama SVD code.
4603  UnitTestCompressedMatrix<Real>();
4604  UnitTestCompressedMatrix2<Real>();
4605  UnitTestExtractCompressedMatrix<Real>();
4606  UnitTestResize<Real>();
4607  UnitTestResizeCopyDataDifferentStrideType<Real>();
4608  UnitTestNonsymmetricPower<Real>();
4609  UnitTestEigSymmetric<Real>();
4610  KALDI_LOG << " Point A";
4611  UnitTestComplexPower<Real>();
4612  UnitTestEig<Real>();
4613  UnitTestEigSp<Real>();
4614  // commenting these out for now-- they test the speed, but take a while.
4615  // UnitTestSplitRadixRealFftSpeed<Real>();
4616  // UnitTestRealFftSpeed<Real>(); // won't exit!/
4617  UnitTestComplexFt<Real>();
4618  KALDI_LOG << " Point B";
4619  UnitTestComplexFft2<Real>();
4620  UnitTestComplexFft<Real>();
4621  UnitTestSplitRadixComplexFft<Real>();
4622  UnitTestSplitRadixComplexFft2<Real>();
4623  UnitTestDct<Real>();
4624  UnitTestRealFft<Real>();
4625  KALDI_LOG << " Point C";
4626  UnitTestSplitRadixRealFft<Real>();
4627  UnitTestSvd<Real>();
4628  UnitTestSvdNodestroy<Real>();
4629  UnitTestSvdJustvec<Real>();
4630  UnitTestSpAddDiagVec<Real, float>();
4631  UnitTestSpAddDiagVec<Real, double>();
4632  UnitTestSpAddVecVec<Real>();
4633  UnitTestSpInvert<Real>();
4634  KALDI_LOG << " Point D";
4635  UnitTestTpInvert<Real>();
4636  UnitTestIo<Real>();
4637  UnitTestIoCross<Real>();
4638  UnitTestHtkIo<Real>();
4639  UnitTestScale<Real>();
4640  UnitTestTrace<Real>();
4641  KALDI_LOG << " Point E";
4642  CholeskyUnitTestTr<Real>();
4643  UnitTestAxpy<Real>();
4644  UnitTestSimple<Real>();
4645  UnitTestMmul<Real>();
4646  UnitTestMmulSym<Real>();
4647  UnitTestVecmul<Real>();
4648  UnitTestInverse<Real>();
4649  UnitTestMulElements<Real>();
4650  UnitTestDotprod<Real>();
4651  // UnitTestSvdVariants<Real>();
4652  UnitTestPower<Real>();
4653  UnitTestPowerAbs<Real>();
4654  UnitTestHeaviside<Real>();
4655  UnitTestCopySp<Real>();
4656  UnitTestDeterminant<Real>();
4657  KALDI_LOG << " Point F";
4658  UnitTestDeterminantSign<Real>();
4659  UnitTestSger<Real>();
4660  UnitTestAddOuterProductPlusMinus<Real>();
4661  UnitTestTraceProduct<Real>();
4662  UnitTestTransposeScatter<Real>();
4663  UnitTestRankNUpdate<Real>();
4664  UnitTestSpVec<Real>();
4665  UnitTestLimitCondInvert<Real>();
4666  KALDI_LOG << " Point G";
4667  UnitTestLimitCond<Real>();
4668  UnitTestMat2Vec<Real>();
4669  UnitTestFloorCeiling<Real>();
4670  KALDI_LOG << " Point H";
4671  UnitTestCopyRowsAndCols<Real>();
4672  UnitTestSpliceRows<Real>();
4673  UnitTestAddSp<Real>();
4674  UnitTestRemoveRow<Real>();
4675  UnitTestRow<Real>();
4676  UnitTestSubvector<Real>();
4677  UnitTestRange<Real>();
4678  UnitTestSimpleForVec<Real>();
4679  UnitTestSetRandn<Real>();
4680  UnitTestSetRandUniform<Real>();
4681  UnitTestVectorMax<Real>();
4682  UnitTestVectorMin<Real>();
4683  UnitTestSimpleForMat<Real>();
4684  UnitTestTanh<Real>();
4685  UnitTestSigmoid<Real>();
4686  UnitTestSoftHinge<Real>();
4687  UnitTestNorm<Real>();
4688  UnitTestCopyCols<Real>();
4689  UnitTestCopyRows<Real>();
4690  UnitTestCopyToRows<Real>();
4691  UnitTestAddRows<Real>();
4692  UnitTestAddToRows<Real>();
4693  UnitTestMul<Real>();
4694  KALDI_LOG << " Point I";
4695  UnitTestSolve<Real>();
4696  UnitTestAddMat2<Real>();
4697  UnitTestSymAddMat2<Real>();
4698  UnitTestAddMatSelf<Real>();
4699  UnitTestMaxMin<Real>();
4700  UnitTestInnerProd<Real>();
4701  UnitTestApplyExpSpecial<Real>();
4702  UnitTestScaleDiag<Real>();
4703  UnitTestSetDiag<Real>();
4704  UnitTestSetRandn<Real>();
4705  KALDI_LOG << " Point J";
4706  UnitTestTraceSpSpLower<Real>();
4707  UnitTestTranspose<Real>();
4708  UnitTestAddVec2Sp<Real>();
4709  UnitTestAddVecToRows<Real>();
4710  UnitTestAddVecToCols<Real>();
4712  UnitTestTp2Sp<Real>();
4713  UnitTestTp2<Real>();
4714  UnitTestAddDiagMat2<Real>();
4715  UnitTestAddDiagMatMat<Real>();
4716  // UnitTestOrthogonalizeRows<Real>();
4717  UnitTestTopEigs<Real>();
4718  UnitTestRandCategorical<Real>();
4719  UnitTestTridiag<Real>();
4720  UnitTestTridiag<Real>();
4721  // SlowMatMul<Real>();
4722  UnitTestAddDiagVecMat<Real>();
4723  UnitTestAddMatDiagVec<Real>();
4724  UnitTestAddMatMatElements<Real>();
4725  UnitTestAddMatMatNans<Real>();
4726  UnitTestAddToDiagMatrix<Real>();
4727  UnitTestAddToDiag<Real>();
4728  UnitTestMaxAbsEig<Real>();
4729  UnitTestMax2<Real>();
4730  UnitTestPca<Real>(full_test);
4731  UnitTestPca2<Real>(full_test);
4732  UnitTestAddVecVec<Real>();
4733  UnitTestReplaceValue<Real>();
4734  // The next one is slow. The upshot is that Eig is up to ten times faster
4735  // than SVD.
4736  // UnitTestSvdSpeed<Real>();
4737  KALDI_LOG << " Point K";
4738  UnitTestTriVecSolver<Real>();
4739 }
4740 
4741 
4742 
4743 }
4744 
4745 
4746 int main() {
4747  using namespace kaldi;
4748  bool full_test = false;
4749  SetVerboseLevel(5);
4750  kaldi::MatrixUnitTest<float>(full_test);
4751  kaldi::MatrixUnitTest<double>(full_test);
4752  KALDI_LOG << "Tests succeeded.";
4753 }
static void SlowMatMul()
void CopyColsFromMat(const MatrixBase< Real > &M)
Performs a column stack of the matrix M.
static void UnitTestAddMatMatElements()
bool ApproxEqual(const VectorBase< Real > &other, float tol=0.01) const
Returns true if ((*this)-other).Norm(2.0) <= tol * (*this).Norm(2.0).
static void UnitTestTraceProduct()
static void UnitTestAddVecVec()
void AddMat2(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const Real beta)
rank-N update: if (transM == kNoTrans) (*this) = beta*(*this) + alpha * M * M^T, or (if transM == kTr...
Definition: sp-matrix.cc:1110
This code computes Goodness of Pronunciation (GOP) and extracts phone-level pronunciation feature for...
Definition: chain.dox:20
static void UnitTestInverse()
double Exp(double x)
Definition: kaldi-math.h:83
static void UnitTestSvdJustvec()
void AddVec2Sp(const Real alpha, const VectorBase< Real > &v, const SpMatrix< Real > &S, const Real beta)
Does *this = beta * *thi + alpha * diag(v) * S * diag(v)
Definition: sp-matrix.cc:922
static void UnitTestHeaviside()
void Add(const Real alpha)
Add a scalar to each element.
static void UnitTestCopyCols()
void DoStep(Real function_value, const VectorBase< Real > &gradient)
The user calls this function to provide the class with the function and gradient info at the point Ge...
bool IsUnit(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:480
void CopyFromVec(const SubVector< OtherReal > &orig)
CopyFromVec just interprets the vector as having the same layout as the packed matrix.
static void UnitTestEig()
Real Min() const
Returns minimum element of matrix.
This class provides a way for switching between double and float types.
Definition: matrix-common.h:84
Real TraceMatMatMat(const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const MatrixBase< Real > &C, MatrixTransposeType transC)
Returns tr(A B C)
bool IsDiagonal(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:465
static void UnitTestMul()
void UnitTestComplexPower()
Packed symetric matrix class.
Definition: matrix-common.h:62
void Write(std::ostream &out, bool binary) const
write to stream.
This class is a wrapper that enables you to store a matrix in one of three forms: either as a Matrix<...
int16 mSampleSize
Sample size.
Definition: kaldi-matrix.h:961
static void UnitTestMaxAbsEig()
bool InitKaldiInputStream(std::istream &is, bool *binary)
Initialize an opened stream for reading by detecting the binary header and.
Definition: io-funcs-inl.h:306
void CopyColFromVec(const VectorBase< Real > &v, const MatrixIndexT col)
Copy vector into specific column of matrix.
void CopyDiagFromVec(const VectorBase< Real > &v)
Copy vector into diagonal of matrix.
void Tanh(const MatrixBase< Real > &src)
Set each element to the tanh of the corresponding element of "src".
bool IsDiagonal(Real cutoff=1.0e-05) const
Returns true if matrix is Diagonal.
MatrixResizeType
Definition: matrix-common.h:37
void CopyRowToVec(MatrixIndexT row, VectorBase< Real > *v) const
Copies row #row of the matrix into vector v.
void CopyColToVec(MatrixIndexT col, VectorBase< Real > *v) const
Copies column #col of the matrix into vector v.
static int32 DoubleFactorial(int32 i)
static void UnitTestAddMatSelf()
This class describes the options for maximizing various quadratic objective functions.
Definition: sp-matrix.h:443
static void CholeskyUnitTestTr()
void Scale(Real c)
double SolveQuadraticProblem(const SpMatrix< double > &H, const VectorBase< double > &g, const SolverOptions &opts, VectorBase< double > *x)
Definition: sp-matrix.cc:635
bool ApproxEqual(const SpMatrix< Real > &other, float tol=0.01) const
Returns true if ((*this)-other).FrobeniusNorm() <= tol*(*this).FrobeniusNorma()
Definition: sp-matrix.cc:525
void Read(std::istream &in, bool binary, bool add=false)
bool Power(Real pow)
The Power method attempts to take the matrix to a power using a method that works in general for frac...
void SymAddMat2(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transA, Real beta)
*this = beta * *this + alpha * M M^T, for symmetric matrices.
static void UnitTestMmulSym()
void DiffTanh(const MatrixBase< Real > &value, const MatrixBase< Real > &diff)
void UnitTestNonsymmetricPower()
void Transpose()
Transpose the matrix.
void Write(std::ostream &out, bool binary) const
static void UnitTestSpInvert()
bool IsZero(Real cutoff=1.0e-06) const
Returns true if matrix is all zeros.
float RandUniform(struct RandomState *state=NULL)
Returns a random number strictly between 0 and 1.
Definition: kaldi-math.h:151
void RealFftInefficient(VectorBase< Real > *v, bool forward)
Inefficient version of Fourier transform, for testing purposes.
static void UnitTestOrthogonalizeRows()
static void UnitTestDeterminantSign()
static void UnitTestDeterminant()
void GetMatrix(Matrix< BaseFloat > *mat) const
Outputs the contents as a matrix.
void Qr(MatrixBase< Real > *Q)
The symmetric QR algorithm.
Definition: qr.cc:409
void CopyDiagFromMat(const MatrixBase< Real > &M)
Extracts the diagonal of the matrix M.
void AddVecVec(const Real alpha, const VectorBase< Real > &v, const VectorBase< Real > &w)
rank-two update, this <– this + alpha (v w&#39; + w v&#39;).
Definition: sp-matrix.cc:1147
void AddMat2Vec(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const VectorBase< Real > &v, const Real beta=0.0)
Extension of rank-N update: this <– beta*this + alpha * M * diag(v) * M^T.
Definition: sp-matrix.cc:1081
MatrixIndexT NumCols() const
Returns number of columns (or zero for empty matrix).
Definition: kaldi-matrix.h:67
Real Cond() const
Returns condition number by computing Svd.
double TraceMat(const MatrixBase< Real > &A)
Returns trace of matrix.
Real SolveQuadraticMatrixProblem(const SpMatrix< Real > &Q, const MatrixBase< Real > &Y, const SpMatrix< Real > &SigmaInv, const SolverOptions &opts, MatrixBase< Real > *M)
Maximizes the auxiliary function : Like a numerically stable version of .
Definition: sp-matrix.cc:729
Base class which provides matrix operations not involving resizing or allocation. ...
Definition: kaldi-matrix.h:49
Real Trace() const
Definition: sp-matrix.cc:171
static void UnitTestRealFftSpeed()
static void UnitTestComplexFt()
static void UnitTestComplexFft2()
void Solve(const TpMatrix< Real > &M, const MatrixTransposeType trans)
If trans == kNoTrans, solves M x = b, where b is the value of *this at input and x is the value of *t...
static void UnitTestMat2Vec()
bool WriteHtk(std::ostream &os, const MatrixBase< Real > &M, HtkHeader htk_hdr)
Real Max() const
Returns maximum element of matrix.
void ApplyCeiling(Real ceil_val, MatrixIndexT *ceiled_count=nullptr)
Applies ceiling to all elements.
Definition: kaldi-vector.h:155
bool ApproxEqual(const MatrixBase< Real > &other, float tol=0.01) const
Returns true if ((*this)-other).FrobeniusNorm() <= tol * (*this).FrobeniusNorm(). ...
void ApplyLogAndCopy(const VectorBase< Real > &v)
Apply natural log to another vector and put result in *this.
static void UnitTestAddVecToRows()
void InvertDouble(Real *logdet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
Definition: sp-matrix.cc:312
int main()
LatticeWeightTpl< FloatType > Plus(const LatticeWeightTpl< FloatType > &w1, const LatticeWeightTpl< FloatType > &w2)
static void UnitTestSimple()
void ComputeDctMatrix(Matrix< Real > *M)
ComputeDctMatrix computes a matrix corresponding to the DCT, such that M * v equals the DCT of vector...
static void InitRandNonsingular(MatrixBase< Real > *M)
void Write(std::ostream &Out, bool binary) const
Writes to C++ stream (option to write in binary).
Real FrobeniusNorm() const
sqrt of sum of square elements.
Definition: sp-matrix.cc:513
static void MatrixUnitTest(bool full_test)
static void UnitTestCompressedMatrix2()
void AddMatSvec(const Real alpha, const MatrixBase< Real > &M, const MatrixTransposeType trans, const VectorBase< Real > &v, const Real beta)
This is as AddMatVec, except optimized for where v contains a lot of zeros.
Real Trace(bool check_square=true) const
Returns trace of matrix.
static void UnitTestRandCategorical()
static void UnitTestRankNUpdate()
static void UnitTestSplitRadixComplexFft()
Real TraceSpMat(const SpMatrix< Real > &A, const MatrixBase< Real > &B)
Returns tr(A B).
Definition: sp-matrix.cc:389
Real * RowData(MatrixIndexT i)
Returns pointer to data for one row (non-const)
Definition: kaldi-matrix.h:87
static void UnitTestInnerProd()
static void UnitTestTridiagonalizeAndQr()
void AddDiagMat2(Real alpha, const MatrixBase< Real > &M, MatrixTransposeType trans=kNoTrans, Real beta=1.0)
Add the diagonal of a matrix times itself: *this = diag(M M^T) + beta * *this (if trans == kNoTrans)...
void AddMat(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transA=kNoTrans)
*this += alpha * M [or M^T]
static void UnitTestPowerAbs()
static void UnitTestEigSp()
float RandGauss(struct RandomState *state=NULL)
Definition: kaldi-math.h:155
kaldi::int32 int32
static void UnitTestSplitRadixRealFft()
void ApplyPow(Real exponent)
Takes matrix to a fraction power via Svd.
Definition: sp-matrix.cc:91
const VectorBase< Real > & GetValue(Real *objf_value=NULL) const
This returns the value of the variable x that has the best objective function so far, and the corresponding objective function value if requested.
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
static void UnitTestHtkIo()
static void UnitTestSvdNodestroy()
A class for storing matrices.
Definition: kaldi-matrix.h:823
void ScaleDiag(const Real alpha)
static void UnitTestMmul()
void SetUnit()
< Set to zero
static void UnitTestSimpleForMat()
void Resize(MatrixIndexT length, MatrixResizeType resize_type=kSetZero)
Set vector to a specified size (can be zero).
void UnitTestIo(bool binary)
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 Write(std::ostream &os, bool binary) const
void CopyFromMat(const MatrixBase< OtherReal > &M, MatrixTransposeType trans=kNoTrans)
Copy given matrix. (no resize is done).
void SwapSparseMatrix(SparseMatrix< BaseFloat > *smat)
Swaps the with the given SparseMatrix.
Real Min() const
Returns the minimum value of any element, or +infinity for the empty vector.
bool IsTridiagonal(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:491
MatrixIndexT NumRows() const
void PlaceNansInGaps(Matrix< Real > *mat)
static void UnitTestRealFft()
static void UnitTestAddVecToCols()
static void UnitTestSpVec()
static void UnitTestResize()
void Read(std::istream &is, bool binary)
static bool approx_equal(Real a, Real b)
static void UnitTestDotprod()
void SetUnit()
Sets to zero, except ones along diagonal [for non-square matrices too].
void Read(std::istream &is, bool binary)
Note: if you write a compressed matrix in text form, it will be read as a regular full matrix...
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.
Definition: sp-matrix.cc:827
void CopyFromSp(const SpMatrix< Real > &other)
Definition: sp-matrix.h:85
Real Norm(Real p) const
Compute the p-th norm of the vector.
void Tridiagonalize(MatrixBase< Real > *Q)
Tridiagonalize the matrix with an orthogonal transformation.
Definition: qr.cc:165
void SetRandUniform()
Sets to numbers uniformly distributed on (0, 1)
void CopyRowFromMat(const MatrixBase< Real > &M, MatrixIndexT row)
Extracts a row of the matrix M.
static void UnitTestMulElements()
static void UnitTestSymAddMat2()
static void UnitTestCopyToRows()
static void UnitTestTriVecSolver()
int32 mSamplePeriod
Sample period.
Definition: kaldi-matrix.h:959
static void UnitTestFloorCeiling()
void SetRandn()
< Set to unit matrix.
void AddSpSp(const Real alpha, const SpMatrix< Real > &A, const SpMatrix< Real > &B, const Real beta)
this <– beta*this + alpha*A*B.
void CopyFromSp(const SpMatrix< OtherReal > &M)
Copy given spmatrix. (no resize is done).
void CopyToMat(MatrixBase< OtherReal > *other, MatrixTransposeType t=kNoTrans) const
Copy to matrix. It must already have the correct size.
void SetVerboseLevel(int32 i)
This should be rarely used, except by programs using Kaldi as library; command-line programs set the ...
Definition: kaldi-error.h:64
static void UnitTestAxpy()
Real LogSumExp(Real prune=-1.0) const
Returns log(sum(exp())) without exp overflow If prune > 0.0, it uses a pruning beam, discarding terms less than (max - prune).
void AddVec2(const Real alpha, const VectorBase< Real > &v)
Add vector : *this = *this + alpha * rv^2 [element-wise squaring].
int ApplyFloor(const SpMatrix< Real > &Floor, Real alpha=1.0, bool verbose=false)
Floors this symmetric matrix to the matrix alpha * Floor, where the matrix Floor is positive definite...
Definition: sp-matrix.cc:536
void ApplyFloor(Real floor_val, MatrixIndexT *floored_count=nullptr)
Applies floor to all elements.
Definition: kaldi-vector.h:149
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...
bool IsZero(Real cutoff=1.0e-05) const
Returns true if matrix is all zeros.
void CopyFromVec(const VectorBase< Real > &v)
Copy data from another vector (must match own size).
void AddToDiag(const Real r)
MatrixIndexT NumCols() const
void SoftHinge(const MatrixBase< Real > &src)
Set each element to y = log(1 + exp(x))
void AddVec2(const Real alpha, const VectorBase< OtherReal > &v)
rank-one update, this <– this + alpha v v&#39;
Definition: sp-matrix.cc:946
void UnitTestAddVecCross()
static void UnitTestTridiag()
void CopyToRows(Real *const *dst) const
Copies row r of this matrix to the array of floats at the location given by dst[r].
void Read(std::istream &in, bool binary, bool add=false)
read from stream.
bool IsSymmetric(Real cutoff=1.0e-05) const
Returns true if matrix is Symmetric.
void Cholesky(const SpMatrix< Real > &orig)
Definition: tp-matrix.cc:88
MatrixIndexT Stride() const
Stride (distance in memory between each row). Will be >= NumCols.
Definition: kaldi-matrix.h:70
static void UnitTestNorm()
int32 MatrixIndexT
Definition: matrix-common.h:98
void Compute(Real *xr, Real *xi, bool forward) const
Definition: srfft.cc:136
static void UnitTestSigmoid()
static void UnitTestAddMatMatNans()
void AddSmat2Sp(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const SpMatrix< Real > &A, const Real beta=0.0)
This is a version of AddMat2Sp specialized for when M is fairly sparse.
Definition: sp-matrix.cc:1026
static void UnitTestReplaceValue()
Real LogPosDefDet() const
Computes log determinant but only for +ve-def matrices (it uses Cholesky).
Definition: sp-matrix.cc:36
Real TraceMatSpMat(const MatrixBase< Real > &A, MatrixTransposeType transA, const SpMatrix< Real > &B, const MatrixBase< Real > &C, MatrixTransposeType transC)
Returns tr(A B C) (A and C may be transposed as specified by transA and transC).
Definition: sp-matrix.cc:415
static void UnitTestEigSymmetric()
static void UnitTestCopyRows()
const SubVector< Real > Row(MatrixIndexT i) const
Return specific row of matrix [const].
Definition: kaldi-matrix.h:188
double Log(double x)
Definition: kaldi-math.h:100
MatrixIndexT LimitCond(Real maxCond=1.0e+5, bool invert=false)
Definition: sp-matrix.cc:606
void Scale(Real alpha)
Multiply each element with a scalar value.
static void UnitTestSplitRadixRealFftSpeed()
void MulElements(const VectorBase< Real > &v)
Multiply element-by-element by another vector.
static void UnitTestExtractCompressedMatrix()
void AddSp(const Real alpha, const SpMatrix< Real > &Ma)
Definition: sp-matrix.h:211
static void UnitTestSplitRadixComplexFft2()
static void UnitTestAddDiagMatMat()
static void UnitTestSvdZero()
void ApplyPowAbs(Real power, bool include_sign=false)
Take the absolute value of all elements of a vector to a power.
void AddVecToRows(const Real alpha, const VectorBase< OtherReal > &v)
[each row of *this] += alpha * v
static void UnitTestComplexFft()
static void UnitTestLimitCond()
static void UnitTestAddSp()
static void UnitTestGeneralMatrix()
struct rnnlm::@11::@12 n
MatrixStrideType
Definition: matrix-common.h:44
static void UnitTestTp2()
static void UnitTestTridiagonalize()
static void UnitTestSvdSpeed()
void RandPosdefSpMatrix(MatrixIndexT dim, SpMatrix< Real > *matrix)
static void UnitTestSetRandn()
void AddMatSp(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const SpMatrix< Real > &B, const Real beta)
this <– beta*this + alpha*A*B.
Definition: kaldi-matrix.h:708
static void UnitTestTanh()
void SetDiag(const Real alpha)
void CopyFromTp(const TpMatrix< OtherReal > &M, MatrixTransposeType trans=kNoTrans)
Copy given tpmatrix. (no resize is done).
void SetRandn()
Sets to random values of a normal distribution.
static void UnitTestPca(bool full_test)
void AddMatMat(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const Real beta)
static void UnitTestAddDiagMat2()
static void UnitTestPca2(bool full_test)
void AddToRows(Real alpha, Real *const *dst) const
For each row r of this matrix, adds it (times alpha) to the array of floats at the location given by ...
void ExtractRowRangeWithPadding(const GeneralMatrix &in, int32 row_offset, int32 num_rows, GeneralMatrix *out)
This function extracts a row-range of a GeneralMatrix and writes as a GeneralMatrix containing the sa...
#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
static void UnitTestAddMat2()
Real Max() const
Returns the maximum value of any element, or -infinity for the empty vector.
void ComputePca(const MatrixBase< Real > &X, MatrixBase< Real > *U, MatrixBase< Real > *A, bool print_eigs, bool exact)
ComputePCA does a PCA computation, using either outer products or inner products, whichever is more e...
double TraceSpSp(const SpMatrix< double > &A, const SpMatrix< double > &B)
Definition: sp-matrix.cc:326
static void UnitTestSpliceRows()
uint16 mSampleKind
Sample kind.
Definition: kaldi-matrix.h:963
static void UnitTestApplyExpSpecial()
void AddMatMatMat(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const MatrixBase< Real > &C, MatrixTransposeType transC, const Real beta)
this <– beta*this + alpha*A*B*C.
Real VecSpVec(const VectorBase< Real > &v1, const SpMatrix< Real > &M, const VectorBase< Real > &v2)
Computes v1^T * M * v2.
Definition: sp-matrix.cc:964
Real TraceSpSpLower(const SpMatrix< Real > &A, const SpMatrix< Real > &B)
Definition: sp-matrix.cc:1171
Packed symetric matrix class.
Definition: matrix-common.h:63
static void UnitTestAddMat2Sp()
void Scale(float alpha)
scales all elements of matrix by alpha.
#define KALDI_WARN
Definition: kaldi-error.h:150
static void UnitTestTrace()
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.
Real * Data()
Returns a pointer to the start of the vector&#39;s data.
Definition: kaldi-vector.h:70
static void UnitTestSubvector()
void AddOuterProductPlusMinus(Real alpha, const VectorBase< Real > &a, const VectorBase< Real > &b, MatrixBase< Real > *plus, MatrixBase< Real > *minus)
static void UnitTestAddDiagVecMat()
static void UnitTestVectorMin()
static void UnitTestIoCross()
static void UnitTestSolve()
static void UnitTestTransposeScatter()
MatrixIndexT Dim() const
Returns the dimension of the vector.
Definition: kaldi-vector.h:64
Real Sum() const
Returns sum of all elements in matrix.
void SetZero()
Sets matrix to zero.
static Real NonOrthogonality(const MatrixBase< Real > &M, MatrixTransposeType transM)
void Scale(Real alpha)
Multiplies all elements by this constant.
static void UnitTestAddMatDiagVec()
static void UnitTestTranspose()
void AddMatVec(const Real alpha, const MatrixBase< Real > &M, const MatrixTransposeType trans, const VectorBase< Real > &v, const Real beta)
Add matrix times vector : this <– beta*this + alpha*M*v.
Definition: kaldi-vector.cc:92
static void UnitTestAddRows()
static void UnitTestDct()
void MulElements(const MatrixBase< Real > &A)
Element by element multiplication with a given matrix.
Real RecentStepLength() const
Returns the average magnitude of the last n steps (but not more than the number we have stored)...
Definition: optimization.cc:55
int Rand(struct RandomState *state)
Definition: kaldi-math.cc:45
Real Sum() const
Returns sum of the elements.
static void UnitTestMax2()
void SetRandn()
Set vector to random normally-distributed noise.
static Real NonUnitness(const SpMatrix< Real > &S)
void AddVecDivVec(Real alpha, const VectorBase< Real > &v, const VectorBase< Real > &r, Real beta)
Add element-by-element quotient of two vectors.
Real TraceMatSpMatSp(const MatrixBase< Real > &A, MatrixTransposeType transA, const SpMatrix< Real > &B, const MatrixBase< Real > &C, MatrixTransposeType transC, const SpMatrix< Real > &D)
Returns tr (A B C D) (A and C may be transposed as specified by transA and transB).
Definition: sp-matrix.cc:438
void Eig(MatrixBase< Real > *P, VectorBase< Real > *eigs_real, VectorBase< Real > *eigs_imag) const
Eigenvalue Decomposition of a square NxN matrix into the form (*this) = P D P^{-1}.
void ComplexFft(VectorBase< Real > *v, bool forward, Vector< Real > *tmp_in)
The function ComplexFft does an Fft on the vector argument v.
void MulColsVec(const VectorBase< Real > &scale)
Equivalent to (*this) = (*this) * diag(scale).
static void UnitTestVectorMax()
void MulRowsVec(const VectorBase< Real > &scale)
Equivalent to (*this) = diag(scale) * (*this).
void AddSmatMat(Real alpha, const SparseMatrix< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, Real beta)
(*this) = alpha * op(A) * B + beta * (*this), where A is sparse.
static void UnitTestVecmul()
void CopyFromMat(const MatrixBase< Real > &orig, SpCopyType copy_type=kTakeMean)
Definition: sp-matrix.cc:112
void CopyFromTp(const TpMatrix< Real > &other)
CopyFromTp copies another triangular matrix into this one.
Definition: tp-matrix.h:104
static void UnitTestAddToDiag()
void CopyColFromMat(const MatrixBase< OtherReal > &M, MatrixIndexT col)
Extracts a column of the matrix M.
void DiffSigmoid(const MatrixBase< Real > &value, const MatrixBase< Real > &diff)
static void UnitTestAddVec2Sp()
static void UnitTestAddToDiagMatrix()
static void UnitTestAddOuterProductPlusMinus()
void AddTpVec(const Real alpha, const TpMatrix< Real > &M, const MatrixTransposeType trans, const VectorBase< Real > &v, const Real beta)
Add triangular matrix times vector: this <– beta*this + alpha*M*v.
A class representing a vector.
Definition: kaldi-vector.h:406
void SetRandn(BaseFloat zero_prob)
Sets up to a pseudo-randomly initialized matrix, with each element zero with probability zero_prob an...
static void UnitTestMaxMin()
static void UnitTestCompressedMatrix()
MatrixIndexT NumRows() const
Returns number of rows (or zero for emtpy matrix).
static void UnitTestScaleDiag()
static void UnitTestTp2Sp()
void AddVecToCols(const Real alpha, const VectorBase< OtherReal > &v)
[each col of *this] += alpha * v
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
int32 g_kaldi_verbose_level
This is set by util/parse-options.
Definition: kaldi-error.cc:46
MatrixIndexT NumRows() const
Returns number of rows (or zero for empty matrix).
Definition: kaldi-matrix.h:64
MatrixIndexT NumRows() const
void Set(Real f)
Set all members of a vector to a specified value.
void AddMat2Sp(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const SpMatrix< Real > &A, const Real beta=0.0)
Extension of rank-N update: this <– beta*this + alpha * M * A * M^T.
Definition: sp-matrix.cc:982
Real FrobeniusNorm() const
Frobenius norm, which is the sqrt of sum of square elements.
Real LogDet(Real *det_sign=NULL) const
Returns logdet of matrix.
void ApplyPow(Real power)
Take all elements of vector to a power.
Definition: kaldi-vector.h:179
static void UnitTestSimpleForVec()
void CreateEigenvalueMatrix(const VectorBase< Real > &re, const VectorBase< Real > &im, MatrixBase< Real > *D)
Creates the eigenvalue matrix D that is part of the decomposition used Matrix::Eig.
static void UnitTestSpAddVecVec()
static void UnitTestTopEigs()
Real ApplySoftMax()
Apply soft-max to the collection of all elements of the matrix and return normalizer (log sum of expo...
void AddSp(const Real alpha, const SpMatrix< OtherReal > &S)
*this += alpha * S
Real Cond() const
Returns maximum ratio of singular values.
Definition: sp-matrix.h:145
void AddVecVec(const Real alpha, const VectorBase< OtherReal > &a, const VectorBase< OtherReal > &b)
*this += alpha * a * b^T
MatrixTransposeType
Definition: matrix-common.h:32
int32 LinearCgd(const LinearCgdOptions &opts, const SpMatrix< Real > &A, const VectorBase< Real > &b, VectorBase< Real > *x)
static void AssertEqual(float a, float b, float relative_tolerance=0.001)
assert abs(a - b) <= relative_tolerance * (abs(a)+abs(b))
Definition: kaldi-math.h:276
#define KALDI_VLOG(v)
Definition: kaldi-error.h:156
static void UnitTestRemoveRow()
void CopyRowFromVec(const VectorBase< Real > &v, const MatrixIndexT row)
Copy vector into specific row of matrix.
static void UnitTestRange()
void CopyRowsFromMat(const MatrixBase< Real > &M)
Performs a row stack of the matrix M.
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
static void UnitTestSvd()
static void UnitTestSetDiag()
void cblas_Xspmv(const float alpha, const int num_rows, const float *Mdata, const float *v, const int v_inc, const float beta, float *y, const int y_inc)
void AddDiagMatMat(Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const MatrixBase< Real > &N, MatrixTransposeType transN, Real beta=1.0)
Add the diagonal of a matrix product: *this = diag(M N), assuming the "trans" arguments are both kNoT...
void Resize(const MatrixIndexT r, const MatrixIndexT c, MatrixResizeType resize_type=kSetZero, MatrixStrideType stride_type=kDefaultStride)
Sets matrix to a specified size (zero is OK as long as both r and c are zero).
Real ApplyLogSoftMax()
Applies log soft-max to vector and returns normalizer (log sum of exponentials).
void RemoveRow(MatrixIndexT i)
Remove a specified row.
static void UnitTestFloorUnit()
void InitKaldiOutputStream(std::ostream &os, bool binary)
InitKaldiOutputStream initializes an opened stream for writing by writing an optional binary header a...
Definition: io-funcs-inl.h:291
void ReplaceValue(Real orig, Real changed)
Set each element to y = (x == orig ? changed : x).
static void UnitTestScale()
This is an implementation of L-BFGS.
Definition: optimization.h:84
static void UnitTestSetRandUniform()
void Svd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt) const
Compute SVD (*this) = U diag(s) Vt.
double Log1p(double x)
Definition: kaldi-math.h:104
void CopyFromMat(const MatrixBase< Real > &M, MatrixTransposeType Trans=kNoTrans)
CopyFromMat copies the lower triangle of M into *this (or the upper triangle, if Trans == kTrans)...
Definition: tp-matrix.cc:117
static void UnitTestLimitCondInvert()
static void UnitTestCopySp()
static void UnitTestSoftHinge()
MatrixIndexT RandCategorical() const
This function returns a random index into this vector, chosen with probability proportional to the co...
static void UnitTestCopyRowsAndCols()
void Resize(MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
Definition: sp-matrix.h:81
int32 mNSamples
Number of samples.
Definition: kaldi-matrix.h:957
void ApplyAbs()
Take absolute value of each of the elements.
void RemoveElement(MatrixIndexT i)
Remove one element and shifts later elements down.
void Invert(Real *log_det=NULL, Real *det_sign=NULL, bool inverse_needed=true)
matrix inverse.
Definition: kaldi-matrix.cc:38
void SymPosSemiDefEig(VectorBase< Real > *s, MatrixBase< Real > *P, Real tolerance=0.001) const
This is the version of SVD that we implement for symmetric positive definite matrices.
Definition: sp-matrix.cc:57
static void UnitTestLinearCgd()
static void UnitTestSger()
Provides a vector abstraction class.
Definition: kaldi-vector.h:41
void AddColSumMat(Real alpha, const MatrixBase< Real > &M, Real beta=1.0)
Does *this = alpha * (sum of columns of M) + beta * *this.
void ApplyFloor(Real floor_val)
Definition: kaldi-matrix.h:354
MatrixIndexT NumCols() const
Returns number of columns (or zero for emtpy matrix).
bool IsUnit(Real cutoff=1.0e-05) const
Returns true if the matrix is all zeros, except for ones on diagonal.
void ApplyPow(Real power)
Definition: kaldi-matrix.h:341
void CopyRowsFromVec(const VectorBase< Real > &v)
This function has two modes of operation.
static void UnitTestAddMatSmat()
void Invert(Real *logdet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
matrix inverse.
Definition: sp-matrix.cc:219
#define KALDI_LOG
Definition: kaldi-error.h:153
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) ...
bool IsZero(Real cutoff=1.0e-05) const
Definition: sp-matrix.cc:507
void CopyToMat(MatrixBase< Real > *mat, MatrixTransposeType trans=kNoTrans) const
Copies contents to matrix.
bool ReadHtk(std::istream &is, Matrix< Real > *M_ptr, HtkHeader *header_ptr)
Extension of the HTK header.
Real MaxAbsEig() const
Returns the maximum of the absolute values of any of the eigenvalues.
Definition: sp-matrix.cc:67
A structure containing the HTK header.
Definition: kaldi-matrix.h:955
Sub-matrix representation.
Definition: kaldi-matrix.h:988
static void UnitTestTpInvert()
void AddDiagVec(const Real alpha, const VectorBase< OtherReal > &v)
diagonal update, this <– this + diag(v)
Definition: sp-matrix.cc:183
void Read(std::istream &in, bool binary, bool add=false)
Read function using C++ streams.
Represents a non-allocating general vector which can be defined as a sub-vector of higher-level vecto...
Definition: kaldi-vector.h:501
void CopyColsFromVec(const VectorBase< Real > &v)
Copies vector into matrix, column-by-column.
const VectorBase< Real > & GetProposedValue() const
This returns the value at which the function wants us to compute the objective function and gradient...
Definition: optimization.h:134
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
static Real NonDiagonalness(const SpMatrix< Real > &S)
Real LogDet(Real *det_sign=NULL) const
Definition: sp-matrix.cc:579
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 ...
static void UnitTestPower()
static void UnitTestSpAddDiagVec()
static void UnitTestLbfgs()
static void UnitTestRow()
void Compute(Real *x, bool forward)
If forward == true, this function transforms from a sequence of N real points to its complex fourier ...
Definition: srfft.cc:356
void SortSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt, bool sort_on_absolute_value)
Function to ensure that SVD is sorted.
static void UnitTestAddToRows()
void Set(Real)
Sets all elements to a specific value.
static void UnitTestFloorChol()
int32 RandInt(int32 min_val, int32 max_val, struct RandomState *state)
Definition: kaldi-math.cc:95
static void UnitTestResizeCopyDataDifferentStrideType()
Make sure that when Resize() is called with resize_type = kCopyData and a stride_type different from ...
Real TraceMatMatMatMat(const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const MatrixBase< Real > &C, MatrixTransposeType transC, const MatrixBase< Real > &D, MatrixTransposeType transD)
Returns tr(A B C D)
static void UnitTestSvdBad()
void Sigmoid(const MatrixBase< Real > &src)
Set each element to the sigmoid of the corresponding element of "src".
void Write(std::ostream &os, bool binary) const
void RealFft(VectorBase< Real > *v, bool forward)
RealFft is a fourier transform of real inputs.
static void UnitTestTraceSpSpLower()
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