sbdsqr_split1

Source: tests/dense/bdsqr.cpp

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

// Containers.
#include <myramath/utility/Number.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/Vector.h>
#include <myramath/dense/DiagonalMatrix.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/DiagonalMatrixRange.h>

// Algorithms.
#include <myramath/dense/diag.h>
#include <myramath/dense/dimm.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/bdsqr.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/threshold.h>

// Reporting.
#include <tests/myratest.h>

using namespace myra;

namespace {

// Implementation detail of all tests.
template<class Precision> void test_detail(Vector<Precision>& B0, Vector<Precision>& B1, Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  // Form B.
  int N = B0.size();
  auto B = Matrix<Precision>::zeros(N,N);
  for (int n = 0; n < N; ++n)
    B(n,n) = B0(n);
  for (int n = 0; n < N-1; ++n)    
    B(n,n+1) = B1(n);
  // Decompose B = U*S*V'
  auto UVt = bdsqr_inplace(B0,B1);
  const Matrix<Precision>& U = UVt.first;
  const Matrix<Precision>& Vt = UVt.second;
  DiagonalMatrix<Precision> S = diag(B0);
  // Check residual error |R|, where R = B-U*S*V'
  Precision R_error = frobenius(B-U*S*Vt)/frobenius(S);
  myra::out() << "  |B-U*S*V'| = " << R_error << std::endl;
  // Check orthogonality of U and V.
  auto I = Matrix<Precision>::identity(N);
  Precision U_error = frobenius(gemm(U,'T',U)-I)/N;
  Precision V_error = frobenius(gemm(Vt,Vt,'T')-I)/N;  
  myra::out() << "  |U'U-I| = " << U_error << std::endl;
  myra::out() << "  |VV'-I| = " << V_error << std::endl;
  // Check the layout assumptions of U (orthogonal columns) and Vt (orthogonal rows).
  int K = N/2;
  auto Uk = U.left(K);
  auto Vk = Vt.top(K);
  auto Ik = Matrix<Precision>::identity(K);  
  Precision Uk_error = frobenius(gemm(Uk,'T',Uk)-Ik)/N;
  Precision Vk_error = frobenius(gemm(Vk,Vk,'T')-Ik)/N;
  myra::out() << "  |Uk'Uk-Ik| = " << Uk_error << std::endl;
  myra::out() << "  |Vk'Vk-Ik| = " << Vk_error << std::endl;
  // Verify singular values are all positive.
  bool positive = true;
  for (int n = 0; n < N; ++n)
    positive &= S(n) >= 0;
  myra::out() << std::boolalpha;
  myra::out() << "  positive S >= 0 = " << positive << std::endl;
  // Verify singular are sorted in descending order.
  bool sorted = true;
  for (int n = 0; n < N-1; ++n)
    sorted &= S(n) >= S(n+1);
  myra::out() << std::boolalpha;
  myra::out() << "  sorted S(n) >= S(n+1) = " << sorted << std::endl;
  myra::out() << std::endl;
  // Apply tests.
  REQUIRE(R_error < tolerance);
  REQUIRE(U_error < tolerance);
  REQUIRE(V_error < tolerance);
  REQUIRE(Uk_error < tolerance);
  REQUIRE(Vk_error < tolerance);
  REQUIRE(positive);
  REQUIRE(sorted);
  }

// Tests a random bidiagonal matrix.
template<class Precision> void test_random(int N, Precision tolerance)
  {
  auto B0 = Vector<Precision>::random(N);
  auto B1 = Vector<Precision>::random(N-1);
  test_detail(B0,B1,tolerance); 
  }

// Tests a bidiagonal matrix with zero diagonal (that is, a zero within B0 vector)
template<class Precision> void test_split0(int N, int S, Precision tolerance)
  {
  auto B0 = Vector<Precision>::random(N);
  auto B1 = Vector<Precision>::random(N-1);
  B0(S) = 0;
  test_detail(B0,B1,tolerance);
  }

// Tests a bidiagonal matrix with zero off-diagonal (that is, a zero within B1 vector)
template<class Precision> void test_split1(int N, int S, Precision tolerance)
  {
  auto B0 = Vector<Precision>::random(N);
  auto B1 = Vector<Precision>::random(N-1);
  B1(S) = 0;
  test_detail(B0,B1,tolerance);
  }

} // namespace

// Random test, float precision.
ADD_TEST("sbdsqr","[dense][lapack]")
  { test_random<float>(50,1.0e-4f); }

// Zero diagonal test, float precision.
ADD_TEST("sbdsqr_split0","[dense][lapack]")
  { test_split0<float>(100,27,1.0e-4f); }

// Zero off-diagonal test, float precision.
ADD_TEST("sbdsqr_split1","[dense][lapack]")
  { test_split1<float>(100,27,1.0e-4f); }

// Stress test for float precision, lots of sizes (N) and lots of trials (T).
ADD_TEST("sbdsqr_stress","[dense][lapack][.]")
  {
  int N = 100;
  int T = 10;
  for (int n = 1; n < N; ++n)
    for (int t = 0; t < T; ++t)
      test_random<float>(n,1.0e-4f);
  }

// Random test, double precision.
ADD_TEST("dbdsqr","[dense][lapack]")
  { test_random<double>(10,1.0e-12); }

// Zero diagonal test, double precision.
ADD_TEST("dbdsqr_split0","[dense][lapack]")
  { test_split0<double>(100,27,1.0e-12); }

// Zero off-diagonal test, double precision.
ADD_TEST("dbdsqr_split1","[dense][lapack]")
  { test_split1<double>(100,27,1.0e-12); }

// Stress test for double precision, lots of sizes (N) and lots of trials (T).
ADD_TEST("dbdsqr_stress","[dense][lapack][.]")
  {
  int N = 100;
  int T = 10;
  for (int n = 1; n < N; ++n)
    for (int t = 0; t < T; ++t)
      test_random<double>(n,1.0e-12);
  }

Results: [PASS]

float
  |B-U*S*V'| = 4.69289e-06
  |U'U-I| = 8.21543e-08
  |VV'-I| = 7.43378e-08
  |Uk'Uk-Ik| = 6.1899e-08
  |Vk'Vk-Ik| = 4.91771e-08
  positive S >= 0 = true
  sorted S(n) >= S(n+1) = true

Go back to Summary of /test programs.