// matrix/matrix-lib-test.cc // Copyright 2009-2012 Microsoft Corporation; Mohit Agarwal; Lukas Burget; // Ondrej Glembek; Saarland University (Author: Arnab Ghoshal); // Go Vivace Inc.; Yanmin Qian; Jan Silovsky; // Johns Hopkins University (Author: Daniel Povey); // Haihua Xu // See ../../COPYING for clarification regarding multiple authors // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. // You may obtain a copy of the License at // http://www.apache.org/licenses/LICENSE-2.0 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE, // MERCHANTABLITY OR NON-INFRINGEMENT. // See the Apache 2 License for the specific language governing permissions and // limitations under the License. #include "matrix/matrix-lib.h" #include #include // This is only needed for UnitTestSvdSpeed, you can // comment it (and that function) out if it causes problems. namespace kaldi { template void RandPosdefSpMatrix(MatrixIndexT dim, SpMatrix *matrix) { MatrixIndexT dim2 = dim + (rand() % 3); // slightly higher-dim. // generate random (non-singular) matrix Matrix tmp(dim, dim2); while (1) { tmp.SetRandn(); if (tmp.Cond() < 100) break; KALDI_LOG << "Condition number of random matrix large " << static_cast(tmp.Cond()) << ", trying again (this is normal)"; } // tmp * tmp^T will give positive definite matrix matrix->AddMat2(1.0, tmp, kNoTrans, 0.0); // Checks that the matrix is indeed pos-def TpMatrix sqrt(dim); sqrt.Cholesky(*matrix); } template static void InitRand(TpMatrix *M) { // start: for (MatrixIndexT i = 0;i < M->NumRows();i++) for (MatrixIndexT j = 0;j<=i;j++) (*M)(i, j) = RandGauss(); } template static void InitRand(Vector *v) { for (MatrixIndexT i = 0;i < v->Dim();i++) (*v)(i) = RandGauss(); } template static void InitRand(VectorBase *v) { for (MatrixIndexT i = 0;i < v->Dim();i++) (*v)(i) = RandGauss(); } template static void InitRand(MatrixBase *M) { start: for (MatrixIndexT i = 0;i < M->NumRows();i++) for (MatrixIndexT j = 0;j < M->NumCols();j++) (*M)(i, j) = RandGauss(); if (M->NumRows() != 0 && M->Cond() > 100) { printf("Condition number of random matrix large %f, trying again (this is normal)\n", (float) M->Cond()); goto start; } } template static void InitRand(SpMatrix *M) { start: for (MatrixIndexT i = 0;i < M->NumRows();i++) for (MatrixIndexT j = 0;j<=i;j++) (*M)(i, j) = RandGauss(); if (M->NumRows() != 0 && M->Cond() > 100) goto start; } /* HERE(2) template static void AssertEqual(const Matrix &A, const Matrix &B, float tol = 0.001) { KALDI_ASSERT(A.NumRows() == B.NumRows()&&A.NumCols() == B.NumCols()); for (MatrixIndexT i = 0;i < A.NumRows();i++) for (MatrixIndexT j = 0;j < A.NumCols();j++) { 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))))); } } template static void AssertEqual(const SpMatrix &A, const SpMatrix &B, float tol = 0.001) { KALDI_ASSERT(A.NumRows() == B.NumRows()&&A.NumCols() == B.NumCols()); for (MatrixIndexT i = 0;i < A.NumRows();i++) for (MatrixIndexT j = 0;j<=i;j++) 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))))); } */ /* HERE: template static bool ApproxEqual(const SpMatrix &A, const SpMatrix &B, Real tol = 0.001) { KALDI_ASSERT(A.NumRows() == B.NumRows()); SpMatrix diff(A); diff.AddSp(1.0, B); Real a = std::max(A.Max(), -A.Min()), b = std::max(B.Max(), -B.Min), d = std::max(diff.Max(), -diff.Min()); return (d <= tol * std::max(a, b)); } */ /* was: template bool ApproxEqual(SpMatrix &A, SpMatrix &B, float tol = 0.001) { KALDI_ASSERT(A.NumRows() == B.NumRows()&&A.NumCols() == B.NumCols()); for (MatrixIndexT i = 0;i < A.NumRows();i++) for (MatrixIndexT j = 0;j<=i;j++) 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; return true; } */ /* HERE template static void AssertEqual(Vector &A, Vector &B, float tol = 0.001) { KALDI_ASSERT(A.Dim() == B.Dim()); for (MatrixIndexT i = 0;i < A.Dim();i++) KALDI_ASSERT(std::abs(A(i)-B(i)) < tol); } template static bool ApproxEqual(Vector &A, Vector &B, float tol = 0.001) { KALDI_ASSERT(A.Dim() == B.Dim()); for (MatrixIndexT i = 0;i < A.Dim();i++) if (std::abs(A(i)-B(i)) > tol) return false; return true; } */ template static void CholeskyUnitTestTr() { for (MatrixIndexT i = 0; i < 5; i++) { MatrixIndexT dimM = 2 + rand() % 10; Matrix M(dimM, dimM); InitRand(&M); SpMatrix S(dimM); S.AddMat2(1.0, M, kNoTrans, 0.0); TpMatrix C(dimM); C.Cholesky(S); Matrix CM(C); TpMatrix Cinv(C); Cinv.Invert(); { Matrix A(C), B(Cinv), AB(dimM, dimM); AB.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0); KALDI_ASSERT(AB.IsUnit()); } SpMatrix S2(dimM); S2.AddMat2(1.0, CM, kNoTrans, 0.0); AssertEqual(S, S2); C.Invert(); Matrix CM2(C); CM2.Invert(); SpMatrix S3(dimM); S3.AddMat2(1.0, CM2, kNoTrans, 0.0); AssertEqual(S, S3); } } template static void SlowMatMul() { MatrixIndexT N = 1000; Matrix M(N, N), P(N, N), Q(N, N); for (MatrixIndexT i = 0; i < 10000; i++) { Q.AddMatMat(1.0, M, kNoTrans, P, kNoTrans, 0.0); } } template static void UnitTestAddToDiagMatrix() { for (int p = 0; p < 2; p++) { MatrixIndexT dimM = 10 + rand() % 2, dimN = 1 + rand() % 5; Matrix M(dimM, dimN), Mcopy(M); BaseFloat alpha = 0.35; M.AddToDiag(alpha); for (MatrixIndexT i = 0; i < dimM && i < dimN; i++) Mcopy(i, i) += alpha; AssertEqual(M, Mcopy); } } template static void UnitTestAddDiagVecMat() { for (int p = 0; p < 2; p++) { MatrixIndexT dimM = 100 + rand() % 255, dimN = 100 + rand() % 255; Real alpha = 0.43243, beta = 1.423; Matrix M(dimM, dimN), N(dimM, dimN); M.SetRandn(); N.SetRandn(); MatrixTransposeType trans = (p % 2 == 0 ? kNoTrans : kTrans); if (trans == kTrans) N.Transpose(); Vector V(dimM); V.SetRandn(); Matrix Mcheck(M); for (int32 r = 0; r < dimM; r++) { SubVector Mcheckrow(Mcheck, r); Vector Nrow(dimN); if (trans == kTrans) Nrow.CopyColFromMat(N, r); else Nrow.CopyFromVec(N.Row(r)); Mcheckrow.Scale(beta); Mcheckrow.AddVec(alpha * V(r), Nrow); } M.AddDiagVecMat(alpha, V, N, trans, beta); AssertEqual(M, Mcheck); KALDI_ASSERT(M.Sum() != 0.0); } } template static void UnitTestAddSp() { for (MatrixIndexT i = 0;i< 10;i++) { MatrixIndexT dimM = 10+rand()%10; SpMatrix S(dimM); InitRand(&S); Matrix M(S), N(S); N.AddSp(2.0, S); M.Scale(3.0); AssertEqual(M, N); } } template static void UnitTestSpAddVec() { for (MatrixIndexT i = 0;i< 10;i++) { BaseFloat alpha = (i<5 ? 1.0 : 0.5); MatrixIndexT dimM = 10+rand()%10; SpMatrix S(dimM); InitRand(&S); SpMatrix T(S); Vector v(dimM); InitRand(&v); S.AddVec(alpha, v); for (MatrixIndexT i = 0; i < dimM; i++) T(i, i) += alpha * v(i); AssertEqual(S, T); } } template static void UnitTestSpAddVecVec() { for (MatrixIndexT i = 0;i< 10;i++) { BaseFloat alpha = (i<5 ? 1.0 : 0.5); MatrixIndexT dimM = 10+rand()%10; SpMatrix S(dimM); InitRand(&S); Matrix T(S); Vector v(dimM), w(dimM); InitRand(&v); InitRand(&w); S.AddVecVec(alpha, v, w); T.AddVecVec(alpha, v, w); T.AddVecVec(alpha, w, v); Matrix U(S); AssertEqual(U, T); } } template static void UnitTestCopyRowsAndCols() { // Test other mode of CopyRowsFromVec, and CopyColsFromVec, // where vector is duplicated. for (MatrixIndexT i = 0; i < 30; i++) { MatrixIndexT dimM = 1 + rand() % 5, dimN = 1 + rand() % 5; Vector w(dimN); // test cross-type version of // CopyRowsFromVec. Vector v(dimM); Matrix M(dimM, dimN), N(dimM, dimN); InitRand(&v); InitRand(&w); M.CopyColsFromVec(v); N.CopyRowsFromVec(w); for (MatrixIndexT r = 0; r < dimM; r++) { for (MatrixIndexT c = 0; c < dimN; c++) { KALDI_ASSERT(M(r, c) == v(r)); KALDI_ASSERT(N(r, c) == w(c)); } } } } template static void UnitTestSpliceRows() { for (MatrixIndexT i = 0;i< 10;i++) { MatrixIndexT dimM = 10+rand()%10, dimN = 10+rand()%10; Vector V(dimM*dimN), V10(dimM*dimN); Vector Vs(std::min(dimM, dimN)), Vs10(std::min(dimM, dimN)); InitRand(&V); Matrix M(dimM, dimN); M.CopyRowsFromVec(V); V10.CopyRowsFromMat(M); AssertEqual(V, V10); for (MatrixIndexT i = 0;i < dimM;i++) for (MatrixIndexT j = 0;j < dimN;j++) KALDI_ASSERT(M(i, j) == V(i*dimN + j)); { Vector V2(dimM), V3(dimM); InitRand(&V2); MatrixIndexT x; M.CopyColFromVec(V2, x = (rand() % dimN)); V3.CopyColFromMat(M, x); AssertEqual(V2, V3); } { Vector V2(dimN), V3(dimN); InitRand(&V2); MatrixIndexT x; M.CopyRowFromVec(V2, x = (rand() % dimM)); V3.CopyRowFromMat(M, x); AssertEqual(V2, V3); } { M.CopyColsFromVec(V); V10.CopyColsFromMat(M); AssertEqual(V, V10); } { M.CopyDiagFromVec(Vs); Vs10.CopyDiagFromMat(M); AssertEqual(Vs, Vs10); } } } template static void UnitTestRemoveRow() { // this is for matrix for (MatrixIndexT p = 0;p< 10;p++) { MatrixIndexT dimM = 10+rand()%10, dimN = 10+rand()%10; Matrix M(dimM, dimN); InitRand(&M); MatrixIndexT i = rand() % dimM; // Row to remove. Matrix N(M); N.RemoveRow(i); for (MatrixIndexT j = 0;j < i;j++) { for (MatrixIndexT k = 0;k < dimN;k++) { KALDI_ASSERT(M(j, k) == N(j, k)); } } for (MatrixIndexT j = i+1;j < dimM;j++) { for (MatrixIndexT k = 0;k < dimN;k++) { KALDI_ASSERT(M(j, k) == N(j-1, k)); } } } // this is for vector for (MatrixIndexT p = 0;p< 10;p++) { MatrixIndexT dimM = 10+rand()%10; Vector V(dimM); InitRand(&V); MatrixIndexT i = rand() % dimM; // Element to remove. Vector N(V); N.RemoveElement(i); for (MatrixIndexT j = 0;j < i;j++) { KALDI_ASSERT(V(j) == N(j)); } for (MatrixIndexT j = i+1;j < dimM;j++) { KALDI_ASSERT(V(j) == N(j-1)); } } } template static void UnitTestSubvector() { Vector V(100); InitRand(&V); Vector V2(100); for (MatrixIndexT i = 0;i < 10;i++) { SubVector tmp(V, i*10, 10); SubVector tmp2(V2, i*10, 10); tmp2.CopyFromVec(tmp); } AssertEqual(V, V2); } // just need this for testing function below. Compute n!! static int32 DoubleFactorial(int32 i) { if (i <= 0) { return 1; } else { return i * DoubleFactorial(i - 2); } } template static void UnitTestSetRandn() { for (MatrixIndexT i = 0; i < 5; i++) { MatrixIndexT rows = 100 + rand() % 50, cols = 100 + rand() % 50; Matrix M(rows, cols); M.SetRandn(); for (MatrixIndexT pow = 1; pow < 5; pow++) { // test moments 1 through 4 of // the distribution. Matrix Mpow(M); Mpow.ApplyPow(pow); Real observed_moment = Mpow.Sum() / (rows * cols); // see http://en.wikipedia.org/wiki/Normal_distribution#Moments, // note that mu = 0 and sigma = 1. Real expected_moment = (pow % 2 == 1 ? 0 : DoubleFactorial(pow - 1)); Real k = 10.0; // This is just a constant we use to give us some wiggle // room before rejecting the distribution... e.g. 10 sigma, // quite approximately. Real allowed_deviation = k * pow / sqrt(static_cast(rows * cols)); // give it a bit more wiggle room for higher powers.. this is quite // unscientific, it would be better to involve the absolute moments or // something like that, and use one of those statistical inequalities, // but it involves the gamma function and it's too much hassle to implement. Real lower_bound = expected_moment - allowed_deviation, upper_bound = expected_moment + allowed_deviation; KALDI_ASSERT(observed_moment >= lower_bound && observed_moment <= upper_bound); } } } template static void UnitTestSetRandUniform() { for (MatrixIndexT i = 0; i < 5; i++) { MatrixIndexT rows = 200 + rand() % 50, cols = 200 + rand() % 50; Matrix M(rows, cols); M.SetRandUniform(); M.Add(-0.5); // we'll be testing the central moments, so // center it around zero first. // Got these moments from http://mathworld.wolfram.com/UniformDistribution.html Vector central_moments(5); central_moments(0) = 0.0; central_moments(1) = 0.0; central_moments(2) = 1.0 / 12; // times (b - a)^2, which equals 1. central_moments(3) = 0.0; central_moments(4) = 1.0 / 80; // times (b - a)^4, which equals 1. for (MatrixIndexT pow = 1; pow < central_moments.Dim(); pow++) { Matrix Mpow(M); Mpow.ApplyPow(pow); Real observed_moment = Mpow.Sum() / (rows * cols); // see http://en.wikipedia.org/wiki/Normal_distribution#Moments, // note that mu = 0 and sigma = 1. Real expected_moment = central_moments(pow); Real k = 10.0; // This is just a constant we use to give us some wiggle // room before rejecting the distribution... e.g. 10 sigma, // quite approximately. Real allowed_deviation = k / sqrt(static_cast(rows * cols)); Real lower_bound = expected_moment - allowed_deviation, upper_bound = expected_moment + allowed_deviation; KALDI_ASSERT(observed_moment >= lower_bound && observed_moment <= upper_bound); } } } template static void UnitTestSimpleForVec() { // testing some simple operaters on vector for (MatrixIndexT i = 0; i < 5; i++) { Vector V(100), V1(100), V2(100), V3(100), V4(100); InitRand(&V); if (i % 2 == 0) { V1.SetZero(); V1.Add(1.0); } else { V1.Set(1.0); } V2.CopyFromVec(V); V3.CopyFromVec(V1); V2.InvertElements(); V3.DivElements(V); V4.CopyFromVec(V3); V4.AddVecDivVec(1.0, V1, V, 0.0); AssertEqual(V3, V2); AssertEqual(V4, V3); V4.MulElements(V); AssertEqual(V4, V1); V3.AddVecVec(1.0, V, V2, 0.0); AssertEqual(V3, V1); Vector V5(V), V6(V1), V7(V1), V8(V); V5.AddVec(1.0, V); V8.Scale(2.0); V6.AddVec2(1.0, V); V7.AddVecVec(1.0, V, V, 1.0); AssertEqual(V6, V7); AssertEqual(V5, V8); } for (MatrixIndexT i = 0; i < 5; i++) { MatrixIndexT dimM = 10 + rand() % 10, dimN = 10 + rand() % 10; Matrix M(dimM, dimN); InitRand(&M); Vector Vr(dimN), Vc(dimM); Vr.AddRowSumMat(0.4, M); Vr.AddRowSumMat(0.3, M, 0.5); // note: 0.3 + 0.4*0.5 = 0.5. Vc.AddColSumMat(0.4, M); Vc.AddColSumMat(0.3, M, 0.5); // note: 0.3 + 0.4*0.5 = 0.5. Vr.Scale(2.0); Vc.Scale(2.0); Vector V2r(dimN), V2c(dimM); for (MatrixIndexT k = 0; k < dimM; k++) { V2r.CopyRowFromMat(M, k); AssertEqual(V2r.Sum(), Vc(k)); } for (MatrixIndexT j = 0; j < dimN; j++) { V2c.CopyColFromMat(M, j); AssertEqual(V2c.Sum(), Vr(j)); } } for (MatrixIndexT i = 0; i < 5; i++) { Vector V(100), V1(100), V2(100); InitRand(&V); V1.CopyFromVec(V); V1.ApplyExp(); Real a = V.LogSumExp(); V2.Set(exp(V.LogSumExp())); V1.DivElements(V2); Real b = V.ApplySoftMax(); AssertEqual(V1, V); AssertEqual(a, b); KALDI_LOG << "a = " << a << ", b = " << b; } for (MatrixIndexT i = 0; i < 5; i++) { MatrixIndexT dimV = 10 + rand() % 10; Real p = 0.5 + RandUniform() * 4.5; Vector V(dimV), V1(dimV), V2(dimV); InitRand(&V); V1.AddVecVec(1.0, V, V, 0.0); // V1:=V.*V. V2.CopyFromVec(V1); AssertEqual(V1.Norm(p), V2.Norm(p)); AssertEqual(sqrt(V1.Sum()), V.Norm(2.0)); } for (MatrixIndexT i = 0; i < 5; i++) { MatrixIndexT dimV = 10 + rand() % 10; Real p = RandUniform() * 1.0e-5; Vector V(dimV); V.Set(p); KALDI_ASSERT(V.IsZero(p)); KALDI_ASSERT(!V.IsZero(p*0.9)); } // Test ApplySoftMax() for matrix. Matrix M(10, 10); InitRand(&M); M.ApplySoftMax(); KALDI_ASSERT( fabs(1.0 - M.Sum()) < 0.01); } template static void UnitTestVectorMax() { int32 dimM = 1 + rand() % 10; Vector V(dimM); V.SetRandn(); Real m = V(0); for (int32 i = 1; i < dimM; i++) m = std::max(m, V(i)); KALDI_ASSERT(m == V.Max()); MatrixIndexT i; KALDI_ASSERT(m == V.Max(&i)); KALDI_ASSERT(m == V(i)); } template static void UnitTestVectorMin() { int32 dimM = 1 + rand() % 10; Vector V(dimM); V.SetRandn(); Real m = V(0); for (int32 i = 1; i < dimM; i++) m = std::min(m, V(i)); KALDI_ASSERT(m == V.Min()); MatrixIndexT i; KALDI_ASSERT(m == V.Min(&i)); KALDI_ASSERT(m == V(i)); } template static void UnitTestReplaceValue(){ MatrixIndexT dim = 10 + rand() % 2; Real orig = 0.1 * (rand() % 100), changed = 0.1 * (rand() % 50); Vector V(dim); V.SetRandn(); V(dim / 2) = orig; Vector V1(V); for (MatrixIndexT i = 0; i < V1.Dim(); i ++) { if (V1(i) == orig) V1(i) = changed; } V.ReplaceValue(orig, changed); AssertEqual(V, V1); } template static void UnitTestNorm() { // test some simple norm properties: scaling. also ApproxEqual test. for (MatrixIndexT p = 0; p < 10; p++) { Real scalar = RandGauss(); if (scalar == 0.0) continue; if (scalar < 0) scalar *= -1.0; MatrixIndexT dimM = 10 + rand() % 10, dimN = 10 + rand() % 10; Matrix M(dimM, dimN); InitRand(&M); SpMatrix S(dimM); InitRand(&S); Vector V(dimN); InitRand(&V); Real Mnorm = M.FrobeniusNorm(), Snorm = S.FrobeniusNorm(), Vnorm1 = V.Norm(1.0), Vnorm2 = V.Norm(2.0), Vnorm3 = V.Norm(3.0); M.Scale(scalar); S.Scale(scalar); V.Scale(scalar); KALDI_ASSERT(ApproxEqual(M.FrobeniusNorm(), Mnorm*scalar)); KALDI_ASSERT(ApproxEqual(S.FrobeniusNorm(), Snorm*scalar)); KALDI_ASSERT(ApproxEqual(V.Norm(1.0), Vnorm1 * scalar)); KALDI_ASSERT(ApproxEqual(V.Norm(2.0), Vnorm2 * scalar)); KALDI_ASSERT(ApproxEqual(V.Norm(3.0), Vnorm3 * scalar)); KALDI_ASSERT(V.ApproxEqual(V)); KALDI_ASSERT(M.ApproxEqual(M)); KALDI_ASSERT(S.ApproxEqual(S)); SpMatrix S2(S); S2.Scale(1.1); KALDI_ASSERT(!S.ApproxEqual(S2)); KALDI_ASSERT(S.ApproxEqual(S2, 0.15)); Matrix M2(M); M2.Scale(1.1); KALDI_ASSERT(!M.ApproxEqual(M2)); KALDI_ASSERT(M.ApproxEqual(M2, 0.15)); Vector V2(V); V2.Scale(1.1); KALDI_ASSERT(!V.ApproxEqual(V2)); KALDI_ASSERT(V.ApproxEqual(V2, 0.15)); } } template static void UnitTestCopyRows() { for (MatrixIndexT p = 0; p < 10; p++) { MatrixIndexT num_rows1 = 10 + rand() % 10, num_rows2 = 10 + rand() % 10, num_cols = 10 + rand() % 10; Matrix M(num_rows1, num_cols); InitRand(&M); Matrix N(num_rows2, num_cols), O(num_rows2, num_cols); std::vector reorder(num_rows2); for (int32 i = 0; i < num_rows2; i++) reorder[i] = -1 + (rand() % (num_rows1 + 1)); N.CopyRows(M, reorder); for (int32 i = 0; i < num_rows2; i++) for (int32 j = 0; j < num_cols; j++) if (reorder[i] < 0) O(i, j) = 0; else O(i, j) = M(reorder[i], j); KALDI_LOG << "M is " << M; KALDI_LOG << "N is " << N; KALDI_LOG << "O is " << O; AssertEqual(N, O); } } template static void UnitTestCopyCols() { for (MatrixIndexT p = 0; p < 10; p++) { MatrixIndexT num_cols1 = 10 + rand() % 10, num_cols2 = 10 + rand() % 10, num_rows = 10 + rand() % 10; Matrix M(num_rows, num_cols1); InitRand(&M); Matrix N(num_rows, num_cols2), O(num_rows, num_cols2); std::vector reorder(num_cols2); for (int32 i = 0; i < num_cols2; i++) reorder[i] = -1 + (rand() % (num_cols1 + 1)); N.CopyCols(M, reorder); for (int32 i = 0; i < num_rows; i++) for (int32 j = 0; j < num_cols2; j++) if (reorder[j] < 0) O(i, j) = 0; else O(i, j) = M(i, reorder[j]); KALDI_LOG << "M is " << M; KALDI_LOG << "N is " << N; KALDI_LOG << "O is " << O; AssertEqual(N, O); } } template static void UnitTestSimpleForMat() { // test some simple operates on all kinds of matrix for (MatrixIndexT p = 0; p < 10; p++) { // for FrobeniousNorm() function MatrixIndexT dimM = 10 + rand() % 10, dimN = 10 + rand() % 10; Matrix M(dimM, dimN); InitRand(&M); { Matrix N(M); Real a = M.LogSumExp(), b = N.ApplySoftMax(); AssertEqual(a, b); AssertEqual(1.0, N.Sum()); } { Matrix N(M); N.Add(2.0); for (MatrixIndexT m = 0; m < dimM; m++) for (MatrixIndexT n = 0; n < dimN; n++) N(m, n) -= 2.0; AssertEqual(M, N); } Matrix N(M), M1(M); M1.MulElements(M); Real tmp1 = sqrt(M1.Sum()); Real tmp2 = N.FrobeniusNorm(); KALDI_ASSERT(std::abs(tmp1 - tmp2) < 0.00001); // for LargestAbsElem() function Vector V(dimM); for (MatrixIndexT i = 0; i < dimM; i++) { for (MatrixIndexT j = 0; j < dimN; j++) { M(i, j) = std::abs(M(i, j)); } std::sort(M.RowData(i), M.RowData(i) + dimN); V(i) = M(i, dimN - 1); } std::sort(V.Data(), V.Data() + dimM); KALDI_ASSERT(std::abs(V(dimM - 1) - N.LargestAbsElem()) < 0.00001); } SpMatrix x(3); x.SetZero(); std::stringstream ss; ss << "DP 3\n"; ss << "4.6863" << '\n'; ss << "3.7062 4.6032" << '\n'; ss << "3.4160 3.7256 5.2474" << '\n'; ss >> x; KALDI_ASSERT(x.IsPosDef() == true); // test IsPosDef() function TpMatrix y(3); y.SetZero(); y.Cholesky(x); KALDI_LOG << "Matrix y is a lower triangular Cholesky decomposition of x:" << '\n'; KALDI_LOG << y << '\n'; // test sp-matrix's LogPosDefDet() function Matrix B(x); Real tmp; Real *DetSign = &tmp; KALDI_ASSERT(std::abs(B.LogDet(DetSign) - x.LogPosDefDet()) < 0.00001); for (MatrixIndexT p = 0; p < 10; p++) { // test for sp and tp matrix's AddSp() and AddTp() function MatrixIndexT dimM = 10 + rand() % 10; SpMatrix S(dimM), S1(dimM); TpMatrix T(dimM), T1(dimM); InitRand(&S); InitRand(&S1); InitRand(&T); InitRand(&T1); Matrix M(S), M1(S1), N(T), N1(T1); S.AddSp(1.0, S1); T.AddTp(1.0, T1); M.AddMat(1.0, M1); N.AddMat(1.0, N1); Matrix S2(S); Matrix T2(T); AssertEqual(S2, M); AssertEqual(T2, N); } for (MatrixIndexT i = 0; i < 10; i++) { // test for sp matrix's AddVec2() function MatrixIndexT dimM = 10 + rand() % 10; SpMatrix M(dimM); Vector V(dimM); InitRand(&M); SpMatrix Md(M); InitRand(&V); SpMatrix Sorig(M); M.AddVec2(0.5, V); Md.AddVec2(static_cast(0.5), V); for (MatrixIndexT i = 0; i < dimM; i++) for (MatrixIndexT j = 0; j < dimM; j++) { KALDI_ASSERT(std::abs(M(i, j) - (Sorig(i, j)+0.5*V(i)*V(j))) < 0.001); KALDI_ASSERT(std::abs(Md(i, j) - (Sorig(i, j)+0.5*V(i)*V(j))) < 0.001); } } } template static void UnitTestRow() { for (MatrixIndexT p = 0;p< 10;p++) { MatrixIndexT dimM = 10+rand()%10, dimN = 10+rand()%10; Matrix M(dimM, dimN); InitRand(&M); MatrixIndexT i = rand() % dimM; // Row to get. Vector V(dimN); V.CopyRowFromMat(M, i); // get row. for (MatrixIndexT k = 0;k < dimN;k++) { AssertEqual(M(i, k), V(k)); } { SpMatrix S(dimN); InitRand(&S); Vector v1(dimN), v2(dimN); Matrix M(S); MatrixIndexT dim2 = rand() % dimN; v1.CopyRowFromSp(S, dim2); v2.CopyRowFromMat(M, dim2); AssertEqual(v1, v2); } MatrixIndexT j = rand() % dimN; // Col to get. Vector W(dimM); W.CopyColFromMat(M, j); // get row. for (MatrixIndexT k = 0;k < dimM;k++) { AssertEqual(M(k, j), W(k)); } } } template static void UnitTestAxpy() { for (MatrixIndexT i = 0;i< 10;i++) { MatrixIndexT dimM = 10+rand()%10, dimN = 10+rand()%10; Matrix M(dimM, dimN), N(dimM, dimN), O(dimN, dimM); InitRand(&M); InitRand(&N); InitRand(&O); Matrix Morig(M); M.AddMat(0.5, N); for (MatrixIndexT i = 0;i < dimM;i++) for (MatrixIndexT j = 0;j < dimN;j++) KALDI_ASSERT(std::abs(M(i, j) - (Morig(i, j)+0.5*N(i, j))) < 0.1); M.CopyFromMat(Morig); M.AddMat(0.5, O, kTrans); for (MatrixIndexT i = 0;i < dimM;i++) for (MatrixIndexT j = 0;j < dimN;j++) KALDI_ASSERT(std::abs(M(i, j) - (Morig(i, j)+0.5*O(j, i))) < 0.1); { float f = 0.5 * (float) (rand() % 3); Matrix N(dimM, dimM); InitRand(&N); Matrix N2(N); Matrix N3(N); N2.AddMat(f, N2, kTrans); N3.AddMat(f, N, kTrans); AssertEqual(N2, N3); // check works same with self as arg. } } } template static void UnitTestCopySp() { // Checking that the various versions of copying // matrix to SpMatrix work the same in the symmetric case. for (MatrixIndexT iter = 0;iter < 5;iter++) { int32 dim = 5 + rand() % 10; SpMatrix S(dim), T(dim); S.SetRandn(); Matrix M(S); T.CopyFromMat(M, kTakeMeanAndCheck); AssertEqual(S, T); T.SetZero(); T.CopyFromMat(M, kTakeMean); AssertEqual(S, T); T.SetZero(); T.CopyFromMat(M, kTakeLower); AssertEqual(S, T); T.SetZero(); T.CopyFromMat(M, kTakeUpper); AssertEqual(S, T); } } template static void UnitTestPower() { for (MatrixIndexT iter = 0;iter < 5;iter++) { // this is for matrix-pow MatrixIndexT dimM = 10 + rand() % 10; Matrix M(dimM, dimM), N(dimM, dimM); InitRand(&M); N.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); // N:=M*M^T. SpMatrix S(dimM); S.CopyFromMat(N); // symmetric so should not crash. S.ApplyPow(0.5); S.ApplyPow(2.0); M.CopyFromSp(S); AssertEqual(M, N); // this is for vector-pow MatrixIndexT dimV = 10 + rand() % 10; Vector V(dimV), V1(dimV), V2(dimV); InitRand(&V); V1.AddVecVec(1.0, V, V, 0.0); // V1:=V.*V. V2.CopyFromVec(V1); V2.ApplyPow(0.5); V2.ApplyPow(2.0); AssertEqual(V1, V2); } } template static void UnitTestHeaviside() { for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixIndexT dimM = 10 + rand() % 10, dimN = 10 + rand() % 10; Matrix M(dimM, dimN), N(dimM, dimN); InitRand(&M); N = M; N.ApplyHeaviside(); for (MatrixIndexT r = 0; r < dimM; r++) { for (MatrixIndexT c = 0; c < dimN; c++) { Real x = M(r, c), y = N(r, c); if (x < 0.0) KALDI_ASSERT(y == 0.0); if (x > 0.0) KALDI_ASSERT(y == 1.0); if (x == 0.0) { KALDI_ASSERT(y >= 0.0 && y <= 1.0); } } } } } template static void UnitTestAddOuterProductPlusMinus() { for (MatrixIndexT iter = 0; iter < 10; iter++) { MatrixIndexT dimM = 10 + rand() % 10; MatrixIndexT dimN = 10 + rand() % 10; Matrix M(dimM, dimN), Plus(dimM, dimN), Minus(dimM, dimN), M2(dimM, dimN); Vector v1(dimM), v2(dimN); for (MatrixIndexT i = 0; i < 5; i++) { InitRand(&v1); InitRand(&v2); Real alpha = 0.333 * ((rand() % 10) - 5); M.AddVecVec(alpha, v1, v2); AddOuterProductPlusMinus(alpha, v1, v2, &Plus, &Minus); M2.SetZero(); M2.AddMat(-1.0, Minus); M2.AddMat(1.0, Plus); AssertEqual(M, M2); KALDI_ASSERT(Minus.Min() >= 0); KALDI_ASSERT(Plus.Min() >= 0); } } } template static void UnitTestSger() { for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixIndexT dimM = 10 + rand() % 10; MatrixIndexT dimN = 10 + rand() % 10; Matrix M(dimM, dimN), M2(dimM, dimN); Vector v1(dimM); InitRand(&v1); Vector v2(dimN); InitRand(&v2); Vector v1d(v1), v2d(v2); M.AddVecVec(1.0f, v1, v2); M2.AddVecVec(1.0f, v1, v2); for (MatrixIndexT m = 0;m < dimM;m++) for (MatrixIndexT n = 0;n < dimN;n++) { KALDI_ASSERT(M(m, n) - v1(m)*v2(n) < 0.01); KALDI_ASSERT(M(m, n) - M2(m, n) < 0.01); } } } template static void UnitTestDeterminant() { // also tests matrix axpy and IsZero() and TraceOfProduct{, T} for (MatrixIndexT iter = 0;iter < 5;iter++) { // First test the 2 det routines are the same int dimM = 10 + rand() % 10; Matrix M(dimM, dimM), N(dimM, dimM); InitRand(&M); N.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); // N:=M*M^T. 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) SpMatrix S(dimM); S.CopyFromMat(N); // symmetric so should not crash. Real logdet = S.LogPosDefDet(); Real logdet2, logdet3, sign2, sign3; logdet2 = N.LogDet(&sign2); logdet3 = S.LogDet(&sign3); KALDI_ASSERT(sign2 == 1.0 && sign3 == 1.0 && std::abs(logdet2-logdet) < 0.1 && std::abs(logdet2 - logdet3) < 0.1); Matrix tmp(dimM, dimM); tmp.SetZero(); tmp.AddMat(1.0, N); tmp.AddMat(-1.0, N, kTrans); // symmetric so tmp should be zero. if ( ! tmp.IsZero(1.0e-04)) { printf("KALDI_ERR: matrix is not zero\n"); KALDI_LOG << tmp; KALDI_ASSERT(0); } Real a = TraceSpSp(S, S), b = TraceMatMat(N, N), c = TraceMatMat(N, N, kTrans); KALDI_ASSERT(std::abs(a-b) < 0.1 && std::abs(b-c) < 0.1); } } template static void UnitTestDeterminantSign() { for (MatrixIndexT iter = 0;iter < 20;iter++) { // First test the 2 det routines are the same int dimM = 10 + rand() % 10; Matrix M(dimM, dimM), N(dimM, dimM); InitRand(&M); N.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); // N:=M*M^T. 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) SpMatrix S(dimM); S.CopyFromMat(N); // symmetric so should not crash. Real logdet = S.LogPosDefDet(); Real logdet2, logdet3, sign2, sign3; logdet2 = N.LogDet(&sign2); logdet3 = S.LogDet(&sign3); KALDI_ASSERT(sign2 == 1.0 && sign3 == 1.0 && std::abs(logdet2-logdet) < 0.01 && std::abs(logdet2 - logdet3) < 0.01); MatrixIndexT num_sign_changes = rand() % 5; for (MatrixIndexT change = 0; change < num_sign_changes; change++) { // Change sign of S's det by flipping one eigenvalue, and N by flipping one row. { Matrix M(S); Matrix U(dimM, dimM), Vt(dimM, dimM); Vector s(dimM); M.Svd(&s, &U, &Vt); // SVD: M = U diag(s) Vt s(rand() % dimM) *= -1; U.MulColsVec(s); M.AddMatMat(1.0, U, kNoTrans, Vt, kNoTrans, 0.0); S.CopyFromMat(M); } // change sign of N: N.Row(rand() % dimM).Scale(-1.0); } // add in a scaling factor too. Real tmp = 1.0 + ((rand() % 5) * 0.01); Real logdet_factor = dimM * log(tmp); N.Scale(tmp); S.Scale(tmp); Real logdet4, logdet5, sign4, sign5; logdet4 = N.LogDet(&sign4); logdet5 = S.LogDet(&sign5); AssertEqual(logdet4, logdet+logdet_factor, 0.01); AssertEqual(logdet5, logdet+logdet_factor, 0.01); if (num_sign_changes % 2 == 0) { KALDI_ASSERT(sign4 == 1); KALDI_ASSERT(sign5 == 1); } else { KALDI_ASSERT(sign4 == -1); KALDI_ASSERT(sign5 == -1); } } } template static void UnitTestSpVec() { // Test conversion back and forth between SpMatrix and Vector. for (MatrixIndexT iter = 0;iter < 1;iter++) { MatrixIndexT dimM =10; // 20 + rand()%10; SpMatrix A(dimM), B(dimM); SubVector vec(A); B.CopyFromVec(vec); AssertEqual(A, B); } } template static void UnitTestSherman() { for (MatrixIndexT iter = 0;iter < 1;iter++) { MatrixIndexT dimM =10; // 20 + rand()%10; MatrixIndexT dimK =2; // 20 + rand()%dimM; Matrix A(dimM, dimM), U(dimM, dimK), V(dimM, dimK); InitRand(&A); InitRand(&U); InitRand(&V); for (MatrixIndexT i = 0;i < (MatrixIndexT)dimM;i++) A(i, i) += 0.0001; // Make sure numerically +ve det-- can by chance be almost singular the way we initialized it (I think) Matrix tmpL(dimM, dimM); tmpL.AddMatMat(1.0, U, kNoTrans, V, kTrans, 0.0); // tmpL =U *V. tmpL.AddMat(1.0, A); tmpL.Invert(); Matrix invA(dimM, dimM); invA.CopyFromMat(A); invA.Invert(); Matrix tt(dimK, dimM), I(dimK, dimK); tt.AddMatMat(1.0, V, kTrans, invA, kNoTrans, 0.0); // tt0 = trans(V) *inv(A) Matrix tt1(dimM, dimK); tt1.AddMatMat(1.0, invA, kNoTrans, U, kNoTrans, 0.0); // tt1=INV A *U Matrix tt2(dimK, dimK); tt2.AddMatMat(1.0, V, kTrans, tt1, kNoTrans, 0.0); for (MatrixIndexT i = 0;i < dimK;i++) for (MatrixIndexT j = 0;j < dimK;j++) { if (i == j) I(i, j) = 1.0; else I(i, j) = 0.0; } tt2.AddMat(1.0, I); // I = identity tt2.Invert(); Matrix tt3(dimK, dimM); tt3.AddMatMat(1.0, tt2, kNoTrans, tt, kNoTrans, 0.0); // tt = tt*tran(V) Matrix tt4(dimM, dimM); tt4.AddMatMat(1.0, U, kNoTrans, tt3, kNoTrans, 0.0); // tt = U*tt Matrix tt5(dimM, dimM); tt5.AddMatMat(1.0, invA, kNoTrans, tt4, kNoTrans, 0.0); // tt = inv(A)*tt Matrix tmpR(dimM, dimM); tmpR.CopyFromMat(invA); tmpR.AddMat(-1.0, tt5); // printf("#################%f###############################.%f##############", tmpL, tmpR); AssertEqual(tmpL, tmpR); // checks whether LHS = RHS or not... } } template static void UnitTestTraceProduct() { for (MatrixIndexT iter = 0;iter < 5;iter++) { // First test the 2 det routines are the same int dimM = 10 + rand() % 10, dimN = 10 + rand() % 10; Matrix M(dimM, dimN), N(dimM, dimN); InitRand(&M); InitRand(&N); Matrix Nt(N, kTrans); Real a = TraceMatMat(M, Nt), b = TraceMatMat(M, N, kTrans); printf("m = %d, n = %d\n", dimM, dimN); KALDI_LOG << a << " " << b << '\n'; KALDI_ASSERT(std::abs(a-b) < 0.1); } } template static void UnitTestSvd() { MatrixIndexT Base = 3, Rand = 2, Iter = 25; for (MatrixIndexT iter = 0;iter < Iter;iter++) { MatrixIndexT dimM = Base + rand() % Rand, dimN = Base + rand() % Rand; Matrix M(dimM, dimN); Matrix U(dimM, std::min(dimM, dimN)), Vt(std::min(dimM, dimN), dimN); Vector s(std::min(dimM, dimN)); InitRand(&M); if (iter < 2) KALDI_LOG << "M " << M; Matrix M2(dimM, dimN); M2.CopyFromMat(M); M.Svd(&s, &U, &Vt); if (iter < 2) { KALDI_LOG << " s " << s; KALDI_LOG << " U " << U; KALDI_LOG << " Vt " << Vt; } MatrixIndexT min_dim = std::min(dimM, dimN); Matrix S(min_dim, min_dim); S.CopyDiagFromVec(s); Matrix Mtmp(dimM, dimN); Mtmp.SetZero(); Mtmp.AddMatMatMat(1.0, U, kNoTrans, S, kNoTrans, Vt, kNoTrans, 0.0); AssertEqual(Mtmp, M2); } } template static void UnitTestSvdBad() { MatrixIndexT N = 20; { Matrix M(N, N); // M.Set(1591.3614306764898); M.Set(1.0); M(0, 0) *= 1.000001; Matrix U(N, N), V(N, N); Vector l(N); M.Svd(&l, &U, &V); } SpMatrix M(N); for(MatrixIndexT i =0; i < N; i++) for(MatrixIndexT j = 0; j <= i; j++) M(i, j) = 1591.3614306764898; M(0, 0) *= 1.00001; M(10, 10) *= 1.00001; Matrix U(N, N); Vector l(N); M.SymPosSemiDefEig(&l, &U); } template static void UnitTestSvdZero() { MatrixIndexT Base = 3, Rand = 2, Iter = 30; for (MatrixIndexT iter = 0;iter < Iter;iter++) { MatrixIndexT dimM = Base + rand() % Rand, dimN = Base + rand() % Rand; // M>=N. Matrix M(dimM, dimN); Matrix U(dimM, dimM), Vt(dimN, dimN); Vector v(std::min(dimM, dimN)); if (iter%2 == 0) M.SetZero(); else M.Unit(); if (iter < 2) KALDI_LOG << "M " << M; Matrix M2(dimM, dimN); M2.CopyFromMat(M); bool ans = M.Svd(&v, &U, &Vt); KALDI_ASSERT(ans); // make sure works with zero matrix. } } template static void UnitTestSvdNodestroy() { MatrixIndexT Base = 3, Rand = 2, Iter = 25; for (MatrixIndexT iter = 0;iter < Iter;iter++) { MatrixIndexT dimN = Base + rand() % Rand, dimM = dimN + rand() % Rand; // M>=N, as required by JAMA Svd. MatrixIndexT minsz = std::min(dimM, dimN); Matrix M(dimM, dimN); Matrix U(dimM, minsz), Vt(minsz, dimN); Vector v(minsz); InitRand(&M); if (iter < 2) KALDI_LOG << "M " << M; M.Svd(&v, &U, &Vt); if (iter < 2) { KALDI_LOG << " v " << v; KALDI_LOG << " U " << U; KALDI_LOG << " Vt " << Vt; } for (MatrixIndexT it = 0;it < 2;it++) { Matrix Mtmp(minsz, minsz); for (MatrixIndexT i = 0;i < v.Dim();i++) Mtmp(i, i) = v(i); Matrix Mtmp2(minsz, dimN); Mtmp2.AddMatMat(1.0, Mtmp, kNoTrans, Vt, kNoTrans, 0.0); Matrix Mtmp3(dimM, dimN); Mtmp3.AddMatMat(1.0, U, kNoTrans, Mtmp2, kNoTrans, 0.0); for (MatrixIndexT i = 0;i < Mtmp.NumRows();i++) { for (MatrixIndexT j = 0;j < Mtmp.NumCols();j++) { KALDI_ASSERT(std::abs(Mtmp3(i, j) - M(i, j)) < 0.0001); } } SortSvd(&v, &U, &Vt); // and re-check... for (MatrixIndexT i = 0; i + 1 < v.Dim(); i++) // check SortSvd is working. KALDI_ASSERT(std::abs(v(i+1)) <= std::abs(v(i))); } } } /* template static void UnitTestSvdVariants() { // just make sure it doesn't crash if we call it but don't want left or right singular vectors. there are KALDI_ASSERTs inside the Svd. #ifndef HAVE_ATLAS MatrixIndexT Base = 10, Rand = 5, Iter = 25; for (MatrixIndexT iter = 0;iter < Iter;iter++) { MatrixIndexT dimM = Base + rand() % Rand, dimN = Base + rand() % Rand; // if (dimM=N. Matrix M(dimM, dimN); Matrix U(dimM, dimM), Vt(dimN, dimN); Vector v(std::min(dimM, dimN)); Matrix Utmp(dimM, 1); Matrix Vttmp(1, dimN); InitRand(&M); M.Svd(v, U, Vttmp, "A", "N"); M.Svd(v, Utmp, Vt, "N", "A"); Matrix U2(dimM, dimM), Vt2(dimN, dimN); Vector v2(std::min(dimM, dimN)); M.Svd(v, U2, Vt2, "A", "A"); AssertEqual(U, U2); AssertEqual(Vt, Vt2); } #endif }*/ template static void UnitTestSvdJustvec() { // Making sure gives same answer if we get just the vector, not the eigs. MatrixIndexT Base = 10, Rand = 5, Iter = 25; for (MatrixIndexT iter = 0;iter < Iter;iter++) { MatrixIndexT dimM = Base + rand() % Rand, dimN = Base + rand() % Rand; // M>=N. MatrixIndexT minsz = std::min(dimM, dimN); Matrix M(dimM, dimN); Matrix U(dimM, minsz), Vt(minsz, dimN); Vector v(minsz); M.Svd(&v, &U, &Vt); Vector v2(minsz); M.Svd(&v2); AssertEqual(v, v2); } } template static void UnitTestEigSymmetric() { for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixIndexT dimM = 20 + rand()%10; SpMatrix S(dimM); InitRand(&S); Matrix M(S); // copy to regular matrix. Matrix P(dimM, dimM); Vector real_eigs(dimM), imag_eigs(dimM); M.Eig(&P, &real_eigs, &imag_eigs); KALDI_ASSERT(imag_eigs.Sum() == 0.0); // should have M = P D P^T Matrix tmp(P); tmp.MulColsVec(real_eigs); // P * eigs Matrix M2(dimM, dimM); M2.AddMatMat(1.0, tmp, kNoTrans, P, kTrans, 0.0); // M2 = tmp * Pinv = P * eigs * P^T AssertEqual(M, M2); // check reconstruction worked. } } template static void UnitTestEig() { for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixIndexT dimM = 1 + iter; /* if (iter < 10) dimM = 1 + rand() % 6; else dimM = 5 + rand()%10; */ Matrix M(dimM, dimM); InitRand(&M); Matrix P(dimM, dimM); Vector real_eigs(dimM), imag_eigs(dimM); M.Eig(&P, &real_eigs, &imag_eigs); { // Check that the eigenvalues match up with the determinant. BaseFloat logdet_check = 0.0, logdet = M.LogDet(); for (MatrixIndexT i = 0; i < dimM ; i++) logdet_check += 0.5 * log(real_eigs(i)*real_eigs(i) + imag_eigs(i)*imag_eigs(i)); AssertEqual(logdet_check, logdet); } Matrix Pinv(P); Pinv.Invert(); Matrix D(dimM, dimM); CreateEigenvalueMatrix(real_eigs, imag_eigs, &D); // check that M = P D P^{-1}. Matrix tmp(dimM, dimM); tmp.AddMatMat(1.0, P, kNoTrans, D, kNoTrans, 0.0); // tmp = P * D Matrix M2(dimM, dimM); M2.AddMatMat(1.0, tmp, kNoTrans, Pinv, kNoTrans, 0.0); // M2 = tmp * Pinv = P * D * Pinv. { // print out some stuff.. Matrix MP(dimM, dimM); MP.AddMatMat(1.0, M, kNoTrans, P, kNoTrans, 0.0); Matrix PD(dimM, dimM); PD.AddMatMat(1.0, P, kNoTrans, D, kNoTrans, 0.0); Matrix PinvMP(dimM, dimM); PinvMP.AddMatMat(1.0, Pinv, kNoTrans, MP, kNoTrans, 0.0); AssertEqual(MP, PD); } AssertEqual(M, M2); // check reconstruction worked. } } template static void UnitTestEigSp() { // Test the Eig function with SpMatrix. // We make sure to test pathological cases, that have // either large zero eigenspaces, or two large // eigenspaces with the same absolute value but +ve // and -ve. Also zero matrix. for (MatrixIndexT iter = 0; iter < 100; iter++) { MatrixIndexT dimM = 1 + (rand() % 10); SpMatrix S(dimM); switch (iter % 5) { case 0: // zero matrix. break; case 1: // general random symmetric matrix. InitRand(&S); break; default: { // start with a random matrix; do decomposition; change the eigenvalues to // try to cover the problematic cases; reconstruct. InitRand(&S); Vector s(dimM); Matrix P(dimM, dimM); S.Eig(&s, &P); // We on purpose create a problematic case where // some eigs are either zero or share a value (+ve or -ve) // with some other eigenvalue. for (MatrixIndexT i = 0; i < dimM; i++) { if (rand() % 10 == 0) s(i) = 0; // set that eig to zero. else if (rand() % 10 < 2) { // set that eig to some other randomly chosen eig, // times random sign. s(i) = (rand()%2 == 0 ? 1 : -1) * s(rand() % dimM); } } // Reconstruct s from the eigenvalues "made problematic." S.AddMat2Vec(1.0, P, kNoTrans, s, 0.0); Real *data = s.Data(); std::sort(data, data+dimM); KALDI_LOG << "Real eigs are: " << s; } } Vector s(dimM); Matrix P(dimM, dimM); S.Eig(&s, &P); KALDI_LOG << "Found eigs are: " << s; SpMatrix S2(dimM); S2.AddMat2Vec(1.0, P, kNoTrans, s, 0.0); { SpMatrix diff(S); diff.AddSp(-1.0, S2); Vector s(dimM); Matrix P(dimM, dimM); diff.Eig(&s, &P); KALDI_LOG << "Eigs of difference are " << s; } KALDI_ASSERT(S.ApproxEqual(S2, 1.0e-03f)); } } // TEMP! template static Real NonOrthogonality(const MatrixBase &M, MatrixTransposeType transM) { SpMatrix S(transM == kTrans ? M.NumCols() : M.NumRows()); S.SetUnit(); S.AddMat2(-1.0, M, transM, 1.0); Real max = 0.0; for (MatrixIndexT i = 0; i < S.NumRows(); i++) for (MatrixIndexT j = 0; j <= i; j++) max = std::max(max, std::abs(S(i, j))); return max; } template static Real NonDiagonalness(const SpMatrix &S) { Real max_diag = 0.0, max_offdiag = 0.0; for (MatrixIndexT i = 0; i < S.NumRows(); i++) for (MatrixIndexT j = 0; j <= i; j++) { if (i == j) { max_diag = std::max(max_diag, std::abs(S(i, j))); } else { max_offdiag = std::max(max_offdiag, std::abs(S(i, j))); } } if (max_diag == 0.0) { if (max_offdiag == 0.0) return 0.0; // perfectly diagonal. else return 1.0; // perfectly non-diagonal. } else { return max_offdiag / max_diag; } } template static Real NonUnitness(const SpMatrix &S) { SpMatrix tmp(S.NumRows()); tmp.SetUnit(); tmp.AddSp(-1.0, S); Real max = 0.0; for (MatrixIndexT i = 0; i < tmp.NumRows(); i++) for (MatrixIndexT j = 0; j <= i; j++) max = std::max(max, std::abs(tmp(i, j))); return max; } template static void UnitTestTridiagonalize() { { float tmp[5]; tmp[4] = 1.0; cblas_Xspmv(1, 0.0, tmp+2, tmp+1, 1, 0.0, tmp+4, 1); KALDI_ASSERT(tmp[4] == 0.0); } for (MatrixIndexT i = 0; i < 4; i++) { MatrixIndexT dim = 40 + rand() % 4; SpMatrix S(dim), S2(dim), R(dim), S3(dim); Matrix Q(dim, dim); InitRand(&S); SpMatrix T(S); T.Tridiagonalize(&Q); KALDI_LOG << "S trace " << S.Trace() << ", T trace " << T.Trace(); //KALDI_LOG << S << "\n" << T; AssertEqual(S.Trace(), T.Trace()); // Also test Trace(). Real ans = 0.0; for (MatrixIndexT j = 0; j < dim; j++) ans += T(j, j); AssertEqual(ans, T.Trace()); AssertEqual(T.LogDet(), S.LogDet()); R.AddMat2(1.0, Q, kNoTrans, 0.0); KALDI_LOG << "Non-unit-ness of R is " << NonUnitness(R); KALDI_ASSERT(R.IsUnit(0.01)); // Check Q is orthogonal. S2.AddMat2Sp(1.0, Q, kTrans, T, 0.0); S3.AddMat2Sp(1.0, Q, kNoTrans, S, 0.0); //KALDI_LOG << "T is " << T; //KALDI_LOG << "S is " << S; //KALDI_LOG << "S2 (should be like S) is " << S2; //KALDI_LOG << "S3 (should be like T) is " << S3; AssertEqual(S, S2); AssertEqual(T, S3); } } template static void UnitTestTridiagonalizeAndQr() { { float tmp[5]; tmp[4] = 1.0; // cblas_sspmv(CblasRowMajor, CblasLower, 1, 0.0, tmp+2, // tmp+1, 1, 0.0, tmp+4, 1); cblas_Xspmv(1, 0.0, tmp+2, tmp+1, 1, 0.0, tmp+4, 1); KALDI_ASSERT(tmp[4] == 0.0); } for (MatrixIndexT i = 0; i < 4; i++) { MatrixIndexT dim = 50 + rand() % 4; SpMatrix S(dim), S2(dim), R(dim), S3(dim), S4(dim); Matrix Q(dim, dim); InitRand(&S); SpMatrix T(S); T.Tridiagonalize(&Q); KALDI_LOG << "S trace " << S.Trace() << ", T trace " << T.Trace(); // KALDI_LOG << S << "\n" << T; AssertEqual(S.Trace(), T.Trace()); // Also test Trace(). Real ans = 0.0; for (MatrixIndexT j = 0; j < dim; j++) ans += T(j, j); AssertEqual(ans, T.Trace()); AssertEqual(T.LogDet(), S.LogDet()); R.AddMat2(1.0, Q, kNoTrans, 0.0); KALDI_LOG << "Non-unit-ness of R after tridiag is " << NonUnitness(R); KALDI_ASSERT(R.IsUnit(0.001)); // Check Q is orthogonal. S2.AddMat2Sp(1.0, Q, kTrans, T, 0.0); S3.AddMat2Sp(1.0, Q, kNoTrans, S, 0.0); //KALDI_LOG << "T is " << T; //KALDI_LOG << "S is " << S; //KALDI_LOG << "S2 (should be like S) is " << S2; //KALDI_LOG << "S3 (should be like T) is " << S3; AssertEqual(S, S2); AssertEqual(T, S3); SpMatrix T2(T); T2.Qr(&Q); R.AddMat2(1.0, Q, kNoTrans, 0.0); KALDI_LOG << "Non-unit-ness of R after QR is " << NonUnitness(R); KALDI_ASSERT(R.IsUnit(0.001)); // Check Q is orthogonal. AssertEqual(T.Trace(), T2.Trace()); KALDI_ASSERT(T2.IsDiagonal()); AssertEqual(T.LogDet(), T2.LogDet()); S4.AddMat2Sp(1.0, Q, kTrans, T2, 0.0); //KALDI_LOG << "S4 (should be like S) is " << S4; AssertEqual(S, S4); } } template static void UnitTestMmul() { for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixIndexT dimM = 20 + rand()%10, dimN = 20 + rand()%10, dimO = 20 + rand()%10; // dims between 10 and 20. // MatrixIndexT dimM = 2, dimN = 3, dimO = 4; Matrix A(dimM, dimN), B(dimN, dimO), C(dimM, dimO); InitRand(&A); InitRand(&B); // // KALDI_LOG <<"a = " << A; // KALDI_LOG<<"B = " << B; C.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0); // C = A * B. // KALDI_LOG << "c = " << C; for (MatrixIndexT i = 0;i < dimM;i++) { for (MatrixIndexT j = 0;j < dimO;j++) { double sum = 0.0; for (MatrixIndexT k = 0;k < dimN;k++) { sum += A(i, k) * B(k, j); } KALDI_ASSERT(std::abs(sum - C(i, j)) < 0.0001); } } } } template static void UnitTestMmulSym() { // Test matrix multiplication on symmetric matrices. for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixIndexT dimM = 20 + rand()%10; Matrix A(dimM, dimM), B(dimM, dimM), C(dimM, dimM), tmp(dimM, dimM), tmp2(dimM, dimM); SpMatrix sA(dimM), sB(dimM), sC(dimM), stmp(dimM); InitRand(&A); InitRand(&B); InitRand(&C); // Make A, B, C symmetric. tmp.CopyFromMat(A); A.AddMat(1.0, tmp, kTrans); tmp.CopyFromMat(B); B.AddMat(1.0, tmp, kTrans); tmp.CopyFromMat(C); C.AddMat(1.0, tmp, kTrans); sA.CopyFromMat(A); sB.CopyFromMat(B); sC.CopyFromMat(C); tmp.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0); // tmp = A * B. tmp2.AddSpSp(1.0, sA, sB, 0.0); // tmp = sA*sB. AssertEqual(tmp, tmp2); tmp2.AddSpSp(1.0, sA, sB, 0.0); // tmp = sA*sB. AssertEqual(tmp, tmp2); } } template static void UnitTestAddVecVec() { for (int32 i = 0; i < 20; i++) { int32 dimM = 5 + rand() % 10, dimN = 5 + rand() % 10; Matrix M(dimM, dimN); M.SetRandn(); Matrix N(M); Vector v(dimM), w(dimN); v.SetRandn(); w.SetRandn(); float alpha = 0.2 * (rand() % 10); M.AddVecVec(alpha, v, w); for (int32 j = 0; j < 20; j++) { int32 dimX = rand() % dimM, dimY = rand() % dimN; AssertEqual(M(dimX, dimY), N(dimX, dimY) + alpha * v(dimX) * w(dimY)); } } } template static void UnitTestVecmul() { for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixTransposeType trans = (iter % 2 == 0 ? kTrans : kNoTrans); MatrixIndexT dimM = 20 + rand()%10, dimN = 20 + rand()%10; // dims between 10 and 20. Real alpha = 0.333, beta = 0.5; Matrix A(dimM, dimN); if (trans == kTrans) A.Transpose(); InitRand(&A); Vector x(dimM), y(dimN); //x.SetRandn(); y.SetRandn(); Vector orig_x(x), x2(x); x.AddMatVec(alpha, A, trans, y, beta); // x = A * y + beta*x. x2.AddMatSvec(alpha, A, trans, y, beta); // x = A * y + beta*x for (MatrixIndexT i = 0; i < dimM; i++) { double sum = beta * orig_x(i); for (MatrixIndexT j = 0; j < dimN; j++) { if (trans == kNoTrans) { sum += alpha * A(i, j) * y(j); } else { sum += alpha * A(j, i) * y(j); } } KALDI_ASSERT(std::abs(sum - x(i)) < 0.0001); } AssertEqual(x, x2); } } template static void UnitTestInverse() { for (MatrixIndexT iter = 0;iter < 10;iter++) { MatrixIndexT dimM = 20 + rand()%10; Matrix A(dimM, dimM), B(dimM, dimM), C(dimM, dimM); InitRand(&A); B.CopyFromMat(A); B.Invert(); C.AddMatMat(1.0, A, kNoTrans, B, kNoTrans, 0.0); // C = A * B. for (MatrixIndexT i = 0;i < dimM;i++) for (MatrixIndexT j = 0;j < dimM;j++) KALDI_ASSERT(std::abs(C(i, j) - (i == j?1.0:0.0)) < 0.1); } } template static void UnitTestMulElements() { for (MatrixIndexT iter = 0; iter < 5; iter++) { MatrixIndexT dimM = 20 + rand()%10, dimN = 20 + rand()%10; Matrix A(dimM, dimN), B(dimM, dimN), C(dimM, dimN); InitRand(&A); InitRand(&B); C.CopyFromMat(A); C.MulElements(B); // C = A .* B (in Matlab, for example). for (MatrixIndexT i = 0;i < dimM;i++) for (MatrixIndexT j = 0;j < dimN;j++) KALDI_ASSERT(std::abs(C(i, j) - (A(i, j)*B(i, j))) < 0.0001); } } template static void UnitTestSpLogExp() { for (MatrixIndexT i = 0; i < 5; i++) { MatrixIndexT dimM = 10 + rand() % 10; Matrix M(dimM, dimM); InitRand(&M); SpMatrix B(dimM); B.AddMat2(1.0, M, kNoTrans, 0.0); // B = M*M^T -> positive definite (since M nonsingular). SpMatrix B2(B); B2.Log(); B2.Exp(); AssertEqual(B, B2); SpMatrix B3(B); B3.Log(); B3.Scale(0.5); B3.Exp(); Matrix sqrt(B3); SpMatrix B4(B.NumRows()); B4.AddMat2(1.0, sqrt, kNoTrans, 0.0); AssertEqual(B, B4); } } template static void UnitTestDotprod() { for (MatrixIndexT iter = 0;iter < 5;iter++) { MatrixIndexT dimM = 200 + rand()%100; Vector v(dimM), w(dimM); InitRand(&v); InitRand(&w); Vector wd(w); Real f = VecVec(w, v), f2 = VecVec(wd, v), f3 = VecVec(v, wd); Real sum = 0.0; for (MatrixIndexT i = 0;i < dimM;i++) sum += v(i)*w(i); KALDI_ASSERT(std::abs(f-sum) < 0.0001); KALDI_ASSERT(std::abs(f2-sum) < 0.0001); KALDI_ASSERT(std::abs(f3-sum) < 0.0001); } } template static void UnitTestResize() { for (size_t i = 0; i < 10; i++) { MatrixIndexT dimM1 = rand() % 10, dimN1 = rand() % 10, dimM2 = rand() % 10, dimN2 = rand() % 10; if (dimM1*dimN1 == 0) dimM1 = dimN1 = 0; if (dimM2*dimN2 == 0) dimM2 = dimN2 = 0; for (MatrixIndexT j = 0; j < 3; j++) { MatrixResizeType resize_type = static_cast(j); Matrix M(dimM1, dimN1); InitRand(&M); Matrix Mcopy(M); Vector v(dimM1); InitRand(&v); Vector vcopy(v); SpMatrix S(dimM1); InitRand(&S); SpMatrix Scopy(S); M.Resize(dimM2, dimN2, resize_type); v.Resize(dimM2, resize_type); S.Resize(dimM2, resize_type); if (resize_type == kSetZero) { KALDI_ASSERT(S.IsZero()); KALDI_ASSERT(v.Sum() == 0.0); KALDI_ASSERT(M.IsZero()); } else if (resize_type == kCopyData) { for (MatrixIndexT i = 0; i < dimM2; i++) { if (i < dimM1) AssertEqual(v(i), vcopy(i)); else KALDI_ASSERT(v(i) == 0); for (MatrixIndexT j = 0; j < dimN2; j++) { if (i < dimM1 && j < dimN1) AssertEqual(M(i, j), Mcopy(i, j)); else AssertEqual(M(i, j), 0.0); } for (MatrixIndexT i2 = 0; i2 < dimM2; i2++) { if (i < dimM1 && i2 < dimM1) AssertEqual(S(i, i2), Scopy(i, i2)); else AssertEqual(S(i, i2), 0.0); } } } } } } template static void UnitTestTp2Sp() { // Tests AddTp2Sp() for (MatrixIndexT iter = 0; iter < 4; iter++) { MatrixIndexT dimM = 10 + rand()%3; TpMatrix T(dimM); InitRand(&T); SpMatrix S(dimM); InitRand(&S); Matrix M(T); for ( MatrixIndexT i = 0; i < dimM; i++) for (MatrixIndexT j = 0; j < dimM; j++) { if (j <= i) AssertEqual(T(i,j), M(i,j)); else AssertEqual(M(i,j), 0.0); } SpMatrix A(dimM), B(dimM); A.AddTp2Sp(0.5, T, (iter < 2 ? kNoTrans : kTrans), S, 0.0); B.AddMat2Sp(0.5, M, (iter < 2 ? kNoTrans : kTrans), S, 0.0); AssertEqual(A, B); } } template static void UnitTestTp2() { // Tests AddTp2() for (MatrixIndexT iter = 0; iter < 4; iter++) { MatrixIndexT dimM = 10 + rand()%3; TpMatrix T(dimM); InitRand(&T); Matrix M(T); SpMatrix A(dimM), B(dimM); A.AddTp2(0.5, T, (iter < 2 ? kNoTrans : kTrans), 0.0); B.AddMat2(0.5, M, (iter < 2 ? kNoTrans : kTrans), 0.0); AssertEqual(A, B); } } template static void UnitTestAddDiagMat2() { for (MatrixIndexT iter = 0; iter < 4; iter++) { MatrixIndexT dimM = 10 + rand() % 3, dimN = 1 + rand() % 4; Vector v(dimM); InitRand(&v); Vector w(v); if (iter % 2 == 1) { Matrix M(dimM, dimN); InitRand(&M); v.AddDiagMat2(0.5, M, kNoTrans, 0.3); Matrix M2(dimM, dimM); M2.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); Vector diag(dimM); diag.CopyDiagFromMat(M2); w.Scale(0.3); w.AddVec(0.5, diag); AssertEqual(w, v); } else { Matrix M(dimN, dimM); InitRand(&M); v.AddDiagMat2(0.5, M, kTrans, 0.3); Matrix M2(dimM, dimM); M2.AddMatMat(1.0, M, kTrans, M, kNoTrans, 0.0); Vector diag(dimM); diag.CopyDiagFromMat(M2); w.Scale(0.3); w.AddVec(0.5, diag); AssertEqual(w, v); } } } template static void UnitTestAddDiagMatMat() { for (MatrixIndexT iter = 0; iter < 4; iter++) { BaseFloat alpha = 0.432 + rand() % 5, beta = 0.043 + rand() % 2; MatrixIndexT dimM = 10 + rand() % 3, dimN = 5 + rand() % 4; Vector v(dimM); Matrix M_orig(dimM, dimN), N_orig(dimN, dimM); M_orig.SetRandn(); N_orig.SetRandn(); MatrixTransposeType transM = (iter % 2 == 0 ? kNoTrans : kTrans); MatrixTransposeType transN = ((iter/2) % 2 == 0 ? kNoTrans : kTrans); Matrix M(M_orig, transM), N(N_orig, transN); InitRand(&v); Vector w(v); w.AddDiagMatMat(alpha, M, transM, N, transN, beta); { Vector w2(v); Matrix MN(dimM, dimM); MN.AddMatMat(1.0, M, transM, N, transN, 0.0); Vector d(dimM); d.CopyDiagFromMat(MN); w2.Scale(beta); w2.AddVec(alpha, d); AssertEqual(w, w2); } } } template static void UnitTestOrthogonalizeRows() { for (MatrixIndexT iter = 0; iter < 100; iter++) { MatrixIndexT dimM = 4 + rand() % 5, dimN = dimM + (rand() % 2); Matrix M(dimM, dimN); for (MatrixIndexT i = 0; i < dimM; i++) { if (rand() % 5 != 0) M.Row(i).SetRandn(); } if (rand() % 2 != 0) { // Multiply by a random square matrix; // keeps it low rank but will be correlated. Harder // test case. Matrix N(dimM, dimM); N.SetRandn(); Matrix tmp(dimM, dimN); tmp.AddMatMat(1.0, N, kNoTrans, M, kNoTrans, 0.0); M.CopyFromMat(tmp); } M.OrthogonalizeRows(); Matrix I(dimM, dimM); I.AddMatMat(1.0, M, kNoTrans, M, kTrans, 0.0); KALDI_ASSERT(I.IsUnit(1.0e-05)); } } template static void UnitTestTransposeScatter() { for (MatrixIndexT iter = 0;iter < 10;iter++) { MatrixIndexT dimA = 10 + rand()%3; MatrixIndexT dimO = 10 + rand()%3; Matrix Af(dimA, dimA); SpMatrix Ap(dimA); Matrix M(dimO, dimA); Matrix Of(dimO, dimO); SpMatrix Op(dimO); Matrix A_MT(dimA, dimO); for (MatrixIndexT i = 0;i < Ap.NumRows();i++) { for (MatrixIndexT j = 0; j<=i; j++) { Ap(i, j) = RandGauss(); } } for (MatrixIndexT i = 0;i < M.NumRows();i++) { for (MatrixIndexT j = 0; j < M.NumCols(); j++) { M(i, j) = RandGauss(); } } /* std::stringstream ss("1 2 3"); ss >> Ap; ss.clear(); ss.str("5 6 7 8 9 10"); ss >> M; */ Af.CopyFromSp(Ap); A_MT.AddMatMat(1.0, Af, kNoTrans, M, kTrans, 0.0); Of.AddMatMat(1.0, M, kNoTrans, A_MT, kNoTrans, 0.0); Op.AddMat2Sp(1.0, M, kNoTrans, Ap, 0.0); // KALDI_LOG << "A" << '\n' << Af << '\n'; // KALDI_LOG << "M" << '\n' << M << '\n'; // KALDI_LOG << "Op" << '\n' << Op << '\n'; for (MatrixIndexT i = 0; i < dimO; i++) { for (MatrixIndexT j = 0; j<=i; j++) { KALDI_ASSERT(std::abs(Of(i, j) - Op(i, j)) < 0.0001); } } A_MT.Resize(dimO, dimA); A_MT.AddMatMat(1.0, Of, kNoTrans, M, kNoTrans, 0.0); Af.AddMatMat(1.0, M, kTrans, A_MT, kNoTrans, 1.0); Ap.AddMat2Sp(1.0, M, kTrans, Op, 1.0); // KALDI_LOG << "Ap" << '\n' << Ap << '\n'; // KALDI_LOG << "Af" << '\n' << Af << '\n'; for (MatrixIndexT i = 0; i < dimA; i++) { for (MatrixIndexT j = 0; j<=i; j++) { KALDI_ASSERT(std::abs(Af(i, j) - Ap(i, j)) < 0.01); } } } } template static void UnitTestRankNUpdate() { for (MatrixIndexT iter = 0;iter < 10;iter++) { MatrixIndexT dimA = 10 + rand()%3; MatrixIndexT dimO = 10 + rand()%3; Matrix Af(dimA, dimA); SpMatrix Ap(dimA); SpMatrix Ap2(dimA); Matrix M(dimO, dimA); InitRand(&M); Matrix N(M, kTrans); Af.AddMatMat(1.0, M, kTrans, M, kNoTrans, 0.0); Ap.AddMat2(1.0, M, kTrans, 0.0); Ap2.AddMat2(1.0, N, kNoTrans, 0.0); Matrix Ap_f(Ap); Matrix Ap2_f(Ap2); AssertEqual(Ap_f, Af); AssertEqual(Ap2_f, Af); } } template static void UnitTestSpInvert() { for (MatrixIndexT i = 0;i < 30;i++) { MatrixIndexT dimM = 6 + rand()%20; SpMatrix M(dimM); for (MatrixIndexT i = 0;i < M.NumRows();i++) for (MatrixIndexT j = 0;j<=i;j++) M(i, j) = RandGauss(); SpMatrix N(dimM); N.CopyFromSp(M); if (rand() % 2 == 0) N.Invert(); else N.InvertDouble(); Matrix Mm(dimM, dimM), Nm(dimM, dimM), Om(dimM, dimM); Mm.CopyFromSp(M); Nm.CopyFromSp(N); Om.AddMatMat(1.0, Mm, kNoTrans, Nm, kNoTrans, 0.0); KALDI_ASSERT(Om.IsUnit( 0.01*dimM )); } } template static void UnitTestTpInvert() { for (MatrixIndexT i = 0;i < 30;i++) { MatrixIndexT dimM = 20 + rand()%10; TpMatrix M(dimM); for (MatrixIndexT i = 0;i < M.NumRows();i++) { for (MatrixIndexT j = 0;j < i;j++) M(i, j) = RandGauss(); M(i, i) = 20 * std::max((Real)0.1, (Real) RandGauss()); // make sure invertible by making it diagonally dominant (-ish) } TpMatrix N(dimM); N.CopyFromTp(M); N.Invert(); TpMatrix O(dimM); Matrix Mm(dimM, dimM), Nm(dimM, dimM), Om(dimM, dimM); Mm.CopyFromTp(M); Nm.CopyFromTp(N); Om.AddMatMat(1.0, Mm, kNoTrans, Nm, kNoTrans, 0.0); KALDI_ASSERT(Om.IsUnit(0.001)); } } template static void UnitTestLimitCondInvert() { for (MatrixIndexT i = 0;i < 10;i++) { MatrixIndexT dimM = 20 + rand()%10; MatrixIndexT dimN = dimM + 1 + rand()%10; SpMatrix B(dimM); Matrix X(dimM, dimN); InitRand(&X); B.AddMat2(1.0, X, kNoTrans, 0.0); // B = X*X^T -> positive definite (almost certainly), since N > M. SpMatrix B2(B); B2.LimitCond(1.0e+10, true); // Will invert. Matrix Bf(B), B2f(B2); Matrix I(dimM, dimM); I.AddMatMat(1.0, Bf, kNoTrans, B2f, kNoTrans, 0.0); KALDI_ASSERT(I.IsUnit(0.1)); } } template static void UnitTestFloorChol() { for (MatrixIndexT i = 0;i < 10;i++) { MatrixIndexT dimM = 20 + rand()%10; MatrixIndexT dimN = 20 + rand()%10; Matrix X(dimM, dimN); InitRand(&X); SpMatrix B(dimM); B.AddMat2(1.0, X, kNoTrans, 0.0); // B = X*X^T -> positive semidefinite. float alpha = (rand() % 10) + 0.5; Matrix M(dimM, dimM); InitRand(&M); SpMatrix C(dimM); C.AddMat2(1.0, M, kNoTrans, 0.0); // C:=M*M^T InitRand(&M); C.AddMat2(1.0, M, kNoTrans, 1.0); // C+=M*M^T (after making new random M) if (i%2 == 0) C.Scale(0.001); // so it's not too much bigger than B (or it's trivial) SpMatrix BFloored(B); BFloored.ApplyFloor(C, alpha); for (MatrixIndexT j = 0;j < 10;j++) { Vector v(dimM); InitRand(&v); Real ip_b = VecSpVec(v, B, v); Real ip_a = VecSpVec(v, BFloored, v); Real ip_c = alpha * VecSpVec(v, C, v); if (i < 3) KALDI_LOG << "alpha = " << alpha << ", ip_a = " << ip_a << " ip_b = " << ip_b << " ip_c = " << ip_c <<'\n'; KALDI_ASSERT(ip_a>=ip_b*0.999 && ip_a>=ip_c*0.999); } } } template static void UnitTestFloorUnit() { for (MatrixIndexT i = 0;i < 5;i++) { MatrixIndexT dimM = 20 + rand()%10; MatrixIndexT dimN = 20 + rand()%10; float floor = (rand() % 10) - 3; Matrix M(dimM, dimN); InitRand(&M); SpMatrix B(dimM); B.AddMat2(1.0, M, kNoTrans, 0.0); // B = M*M^T -> positive semidefinite. SpMatrix BFloored(B); BFloored.ApplyFloor(floor); Vector s(dimM); Matrix P(dimM, dimM); B.SymPosSemiDefEig(&s, &P); Vector