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.