geqpf_rank

Source: tests/dense/geqpf.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/Vector.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/intRange.h>
#include <vector>

// Algorithms.
#include <myramath/dense/triu.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/geqpf.h>
#include <myramath/dense/orgqr.h>
#include <myramath/dense/swaps.h>
#include <myramath/dense/frobenius.h>

// Reporting.
#include <tests/myratest.h>
#include <myramath/utility/stlprint.h>

using namespace myra;
using namespace myra_stlprint;

namespace {

// Test for random tall/skinny A, make sure R is sorted.
template<class Number> void test1(typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  int I = 40;
  int J = 20;
  auto A = Matrix<Number>::random(I,J);
  // Factor P*A = Q*R inplace.
  Matrix<Number> QR = A;
  auto P_tau = geqpf_inplace(QR);
  auto P = P_tau.first;
  auto tau = P_tau.second;
  // Extract R and form Q explicitly. 
  Matrix<Number> R = triu(QR);
  MatrixRange<Number> Q = orgqr_inplace(QR,tau);
  // Verify A = P'*Q*R
  Matrix<Number> B = Q*R;
  swap_columns(perm2swaps(P),B);
  Precision error = frobenius(A-B) / frobenius(A);
  myra::out() << "|A-P'*Q*R| = " << error << std::endl;
  REQUIRE(error < tolerance);
  // Verify the diagonal entries of R are sorted by modulus (big to small).
  bool sorted = true;
  int K = std::min(I,J);
  for (int k = 1; k < K; ++k)
    if ( scalar_norm1(R(k,k)) > scalar_norm1(R(k-1,k-1)) )
      sorted = false;
  REQUIRE(sorted);
  }

// Test for random short/fat A, make sure R is sorted.
template<class Number> void test2(typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  int I = 30;
  int J = 50;
  auto A = Matrix<Number>::random(I,J);
  // Factor P*A = Q*R inplace.
  Matrix<Number> QR = A;
  auto P_tau = geqpf_inplace(QR);
  auto P = P_tau.first;
  auto tau = P_tau.second;
  // Extract R and form Q explicitly. 
  Matrix<Number> R = triu(QR);
  MatrixRange<Number> Q = orgqr_inplace(QR,tau);
  // Verify A = P'*Q*R
  Matrix<Number> B = Q*R;
  swap_columns(perm2swaps(P),B);
  Precision error = frobenius(A-B) / frobenius(A);
  myra::out() << "|A-P'*Q*R| = " << error << std::endl;
  REQUIRE(error < tolerance);
  // Verify the diagonal entries of R are sorted by modulus (big to small).
  bool sorted = true;
  int K = std::min(I,J);
  for (int k = 1; k < K; ++k)
    if ( scalar_norm1(R(k,k)) > scalar_norm1(R(k-1,k-1)) )
      sorted = false;
  REQUIRE(sorted);
  }

// Test for rank-deficient tall/skinny A, make sure its rank is revealed.
template<class Number> void test3(typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  int I = 40;
  int J = 10;
  auto C = Matrix<Number>::random(I,J);
  auto M = Matrix<Number>::random(J,2*J);
  auto A = gemm(C,M);
   // Factor P*A = Q*R inplace. Note A is size Ix2J, but only has J independent columns.
  Matrix<Number> Q = A;
  auto P_tau = geqpf_inplace(Q);
  auto P = P_tau.first;
  auto tau = P_tau.second;
  // Extract R and form Q explicitly. 
  Matrix<Number> R = triu(Q);
  orgqr_inplace(Q,tau);
  // Verify A = P'*Q*R
  Matrix<Number> B = Q*R;
  swap_columns(perm2swaps(P),B);
  Precision error = frobenius(A-B) / frobenius(A);
  myra::out() << "|A-P'*Q*R| = " << error << std::endl;
  REQUIRE(error < tolerance);
  // Verify the rank of A is revealed, should be large gap between R(J-1,J-1) and R(J,J)
  Number r00 = R(J-1,J-1);
  Number r11 = R(J,J);
  myra::out() << "|r00| = " << std::abs(r00) << std::endl; 
  myra::out() << "|r11| = " << std::abs(r11) << std::endl;
  REQUIRE( std::abs(r00)*tolerance > std::abs(r11) );
  // Check norms of R00 and R11, should be a large gap between them, too.
  auto R00 = triu( R.top(J).left(J) );
  auto R11 = triu( R.bottom(J).right(J) );
  myra::out() << "|R00| = " << frobenius(R00) << std::endl; 
  myra::out() << "|R11| = " << frobenius(R11) << std::endl;
  REQUIRE( frobenius(R00)*tolerance > frobenius(R11) );  
  }

} // namespace

ADD_TEST("geqpf_ts_order","[dense][lapack]")  
  {
  test1<NumberS>(1.0e-4f);
  test1<NumberD>(1.0e-8);
  test1<NumberC>(1.0e-4f);
  test1<NumberZ>(1.0e-8);
  }

ADD_TEST("geqpf_sf_order","[dense][lapack]")  
  {
  test2<NumberS>(1.0e-4f);
  test2<NumberD>(1.0e-8);
  test2<NumberC>(1.0e-4f);
  test2<NumberZ>(1.0e-8);
  }

ADD_TEST("geqpf_rank","[dense][lapack]")  
  {
  test3<NumberS>(1.0e-4f);
  test3<NumberD>(1.0e-8);
  test3<NumberC>(1.0e-4f);
  test3<NumberZ>(1.0e-8);
  }

Results: [PASS]

|A-P'*Q*R| = 1.93144e-07
|r00| = 1.73843
|r11| = 1.97701e-06
|R00| = 24.1719
|R11| = 4.01636e-06
|A-P'*Q*R| = 4.49762e-16
|r00| = 1.89254
|r11| = 3.66976e-15
|R00| = 24.0342
|R11| = 9.08485e-15
|A-P'*Q*R| = 2.29179e-07
|r00| = 3.98357
|r11| = 3.58438e-06
|R00| = 43.324
|R11| = 8.74602e-06
|A-P'*Q*R| = 3.4308e-16
|r00| = 4.158
|r11| = 5.79689e-15
|R00| = 43.1554
|R11| = 1.47554e-14

Go back to Summary of /test programs.