test_svd.cc

#include "libmesh/libmesh.h"
#include "libmesh/dense_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/getpot.h"
#include "libmesh/perf_log.h"

#include <Eigen/SVD>

using namespace libMesh;

// A handy typedef for the type of Matrix used.  For single-precision
// calculations, use Eigen::MatrixXf.
// typedef Eigen::MatrixXf EigenMatrixType;
typedef Eigen::MatrixXd EigenMatrixType;
typedef Eigen::VectorXd EigenVectorType;

// A "container" for the parts of a libMesh SVD.  Makes it easier to
// pass around like the Eigen one.
class SVDContainer
{
public:
  DenseMatrix<Number> U, VT;
  DenseVector<Number> sigma;
};

// A templated helper function can fill either a DenseMatrix<T> or an
// EigenMatrixType via duck typing.  The test_num controls which test
// matrix is actually created.  Takes a reference to a GetPot command_line
// object so the matrix can be filled according to command line options.
template <typename T>
void fill_matrix(T & matrix, GetPot & command_line);

// Computes (and logs the performance of) computing the SVD
// decomposition of A.  No least-squares solve is performed.
SVDContainer libmesh_compute_svd(DenseMatrix<Number> & A);

// Compute the SVD of A_Eigen and handle performance logging.
Eigen::JacobiSVD<EigenMatrixType> eigen_compute_svd(const EigenMatrixType & A_Eigen);

// Compare SVDs computed by libmesh and Eigen.  If the singular values
// are unique, they should be in the same (descending) order.
void compare_svd(const SVDContainer & libmesh_svd,
                 const Eigen::JacobiSVD<EigenMatrixType> & svd);



////////////////////////////////////////////////////////////////////////////////
// Main program
////////////////////////////////////////////////////////////////////////////////
int main (int argc, char** argv)
{
  LibMeshInit init (argc, argv);

  // Create a GetPot object to parse the command line
  GetPot command_line (argc, argv);

  // libMesh SVD tests
  DenseMatrix<Number> A;
  fill_matrix(A, command_line);
  SVDContainer libmesh_svd = libmesh_compute_svd(A);

  // Eigen SVD tests
  EigenMatrixType A_Eigen;
  fill_matrix(A_Eigen, command_line);
  Eigen::JacobiSVD<EigenMatrixType> eigen_svd = eigen_compute_svd(A_Eigen);

  // Check differences between libmesh and Eigen SVDs.
  compare_svd(libmesh_svd, eigen_svd);

  return 0;
}



template <typename T>
void fill_matrix(T & matrix, GetPot & command_line)
{
  // Read the test number from the command line.
  unsigned test_num = 1;
  if (command_line.search(1, "--test-num"))
    test_num = command_line.next(test_num);

  switch (test_num)
    {
    case 1:
      matrix.resize(3, 2);
      matrix(0,0) = 2.0; matrix(0,1) = 0.0;
      matrix(1,0) = 0.0; matrix(1,1) = 1.0;
      matrix(2,0) = 0.0; matrix(2,1) = 0.0;
      break;

      // Eigen seems to get a particularly innacurate answer in this case...
    case 2:
      matrix.resize(3, 2);
      matrix(0,0) = 1.0; matrix(0,1) = 2.0;
      matrix(1,0) = 3.0; matrix(1,1) = 4.0;
      matrix(2,0) = 5.0; matrix(2,1) = 6.0;
      break;

    case 3:
      matrix.resize(2, 3);
      matrix(0,0) = 1.0; matrix(0,1) = 2.0; matrix(0,2) = 3.0;
      matrix(1,0) = 4.0; matrix(1,1) = 5.0; matrix(1,2) = 6.0;
      break;

    case 4:
      {
        // You can control the size of the problem by changing nrows and ncols.
        unsigned nrows=50, ncols=6;

        if (command_line.search(1, "--nrows"))
          nrows = command_line.next(nrows);
        if (command_line.search(1, "--ncols"))
          ncols = command_line.next(ncols);

        matrix.resize(nrows, ncols);

        unsigned ctr=1;
        for (unsigned i=0; i<nrows; ++i)
          for (unsigned j=0; j<ncols; ++j)
            matrix(i,j) = static_cast<Real>(ctr++);
        break;
      }
    default:
      libmesh_error_msg("Invalid test_num " << test_num << " requested.");
    }
}



SVDContainer libmesh_compute_svd(DenseMatrix<Number> & A)
{
  PerfLog perf_log("libMesh Build SVD");

  // Note: A.svd() destroys the contents of A, so let's keep a copy!
  DenseMatrix<Number> A_copy = A;

  // Eventual return value
  SVDContainer svd;

  // Compute the SVD
  perf_log.push("libMesh SVD");
  A.svd(svd.sigma, svd.U, svd.VT);
  perf_log.pop("libMesh SVD");

  // Print the results if smaller than 10x10
  if (A_copy.m() <= 10)
    {
      libMesh::out << "## Libmesh DenseMatrix SVD" << std::endl;
      libMesh::out << "A = " << std::endl << A_copy << std::endl;
      libMesh::out << "U = " << std::endl;
      svd.U.print_scientific(libMesh::out, /*precision=*/15);
      libMesh::out << "VT = " << std::endl;
      svd.VT.print_scientific(libMesh::out, /*precision=*/15);
      libMesh::out << "sigma = " << std::endl;
      svd.sigma.print_scientific(libMesh::out, /*precision=*/15);
    }
  
  // Verify results
  DenseMatrix<Number> sigma_matrix(svd.sigma.size(), svd.sigma.size());
  for (unsigned i=0; i<svd.sigma.size(); ++i)
    sigma_matrix(i,i) = svd.sigma(i);

  DenseMatrix<Number> A_computed = svd.U;
  A_computed.right_multiply(sigma_matrix);
  A_computed.right_multiply(svd.VT);

  DenseMatrix<Number> A_diff = A_copy;
  A_diff -= A_computed;
  Real linfty_norm = A_diff.linfty_norm();
  libMesh::out << "libMesh: A_diff.linfty_norm()=" << linfty_norm << std::endl;

  return svd;
}






Eigen::JacobiSVD<EigenMatrixType> eigen_compute_svd(const EigenMatrixType & A_Eigen)
{
  PerfLog perf_log("Eigen Build SVD");

  // Compute the SVD.  There is a thread about how slow the SVD in Eigen is here:
  // https://forum.kde.org/viewtopic.php?f=74&t=102088
  perf_log.push("Eigen SVD");
  Eigen::JacobiSVD<EigenMatrixType> svd(A_Eigen, Eigen::ComputeThinU | Eigen::ComputeThinV);
  perf_log.pop("Eigen SVD");

  // Print the results if smaller than 10x10
  if (A_Eigen.rows() <= 10)
    {
      libMesh::out << "\n## Eigen SVD" << std::endl;
      libMesh::out << "A = " << std::endl << A_Eigen << std::endl << std::endl;
      libMesh::out << "U = " << std::endl << svd.matrixU() << std::endl << std::endl;
      libMesh::out << "VT = " << std::endl << svd.matrixV().transpose() << std::endl << std::endl;
      libMesh::out << "sigma = " << std::endl << svd.singularValues() << std::endl;
    }
  
  // Verify the results
  unsigned n = svd.singularValues().size();
  EigenMatrixType sigma_matrix(n, n);
  for (unsigned i=0; i<n; ++i)
    sigma_matrix(i,i) = svd.singularValues()(i);

  // Compute A using the U, S, V components.
  EigenMatrixType A_computed = svd.matrixU() * sigma_matrix * svd.matrixV().transpose();

  // Compute the norm of the difference ||A - A_computed||.
  Real linfty_norm = (A_Eigen-A_computed).lpNorm<Eigen::Infinity>();
  libMesh::out << "Eigen: A_diff.lpNorm<Infinity>()=" << linfty_norm << std::endl;

  return svd;
}



// A macro to check a user-defined condition and error if failure or
// print a success message on success.
#define check_condition(condition, success_msg, fail_msg) \
  do {                                                    \
    if (!(condition))                                     \
      libmesh_error_msg(fail_msg);                        \
    else                                                  \
      libMesh::out << success_msg << std::endl;           \
  } while(0)



void compare_svd(const SVDContainer & libmesh_svd,
                 const Eigen::JacobiSVD<EigenMatrixType> & svd)
{
  libMesh::out << "\nComparing Results:" << std::endl;

  check_condition(libmesh_svd.U.m() == svd.matrixU().rows(), "Number of rows of U match.", "Number of rows of U do not match.");
  check_condition(libmesh_svd.U.n() == svd.matrixU().cols(), "Number of cols of U match.", "Number of cols of U do not match.");
  check_condition(libmesh_svd.VT.m() == svd.matrixV().transpose().rows(), "Number of rows of VT match.", "Number of rows of VT do not match.");
  check_condition(libmesh_svd.VT.n() == svd.matrixV().transpose().cols(), "Number of cols of VT match.", "Number of cols of VT do not match.");
  check_condition(libmesh_svd.sigma.size() == svd.singularValues().size(), "Number of singualar values match.", "Number of singualar values do not match.");

  // Compare the entries of U and V between the two cases.
  //
  // I don't think it makes sense to compare the U and V vectors,
  // especially when there are very small singualar values, as they
  // may be in completely different order...
  //  {
  //    Real
  //      abs_diff = 0.,
  //      abs_sum = 0.;
  //
  //    for (unsigned i=0; i<libmesh_svd.U.m(); ++i)
  //      for (unsigned j=0; j<libmesh_svd.U.n(); ++j)
  //        {
  //          abs_diff += std::abs(libmesh_svd.U(i,j) - svd.matrixU()(i,j));
  //          abs_sum += std::abs(libmesh_svd.U(i,j) + svd.matrixU()(i,j));
  //        }
  //
  //    Real min_diff = std::min(abs_diff, abs_sum);
  //    libMesh::out << "Norm of difference in U entries: " << min_diff << std::endl;
  //  }
  //  {
  //    Real
  //      abs_diff = 0.,
  //      abs_sum = 0.;
  //
  //    for (unsigned i=0; i<libmesh_svd.VT.m(); ++i)
  //      for (unsigned j=0; j<libmesh_svd.VT.n(); ++j)
  //        {
  //          abs_diff += std::abs(libmesh_svd.VT(i,j) - svd.matrixV()(j,i));
  //          abs_sum += std::abs(libmesh_svd.VT(i,j) + svd.matrixV()(j,i));
  //        }
  //
  //    Real min_diff = std::min(abs_diff, abs_sum);
  //    libMesh::out << "Norm of difference in V^T entries: " << min_diff << std::endl;
  //  }
  {
    Real diff = 0.;
    for (unsigned i=0; i<libmesh_svd.sigma.size(); ++i)
      diff += std::abs(libmesh_svd.sigma(i) - svd.singularValues()(i));
    libMesh::out << "Norm of difference in Sigma entries: " << diff << std::endl;
  }
}