test_svd_least_squares.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 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);

// Verify that the libmesh svd_solve() routine is working.
void libmesh_solve_least_squares(const DenseMatrix<Number> & A);

// Solve a least-squares problem using Eigen.  Requires computation of the SVD.
void eigen_solve_least_squares(const EigenMatrixType & A_Eigen);




////////////////////////////////////////////////////////////////////////////////
// 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);
  libmesh_solve_least_squares(A);

  // Eigen SVD tests
  EigenMatrixType A_Eigen;
  fill_matrix(A_Eigen, command_line);
  eigen_solve_least_squares(A_Eigen);

  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.");
    }
}






void libmesh_solve_least_squares(const DenseMatrix<Number> & A)
{
  PerfLog perf_log("libMesh Least-Squares solve");

  // Construct a rhs vector by multiplying the system matrix with a
  // vector of all 1's.  Then perform an SVD solve using this vector
  // as a rhs, and compare the result to the known solution.
  //
  // Note: if the original matrix has more rows than columns this makes sense
  // (more equations than unknowns).  If the matrix has more columns
  // than rows, there is most likely no solution.
  if (A.m() >= A.n())
    {
      perf_log.push("build soln");
      DenseVector<Number> soln(A.n());
      for (unsigned i=0; i<soln.size(); ++i)
        soln(i) = 1.;
      perf_log.pop("build soln");

      // Compute rhs
      perf_log.push("build rhs");
      DenseVector<Number> rhs;
      A.vector_mult(rhs, soln); // rhs = A * soln;
      perf_log.pop("build rhs");

      // Solve the system of Eqns A*x=rhs for x in the least-squares sense.  Note: A is
      // *not* destroyed by svd_solve(), this function could probably be marked const.
      // For the test matrices that have increasing integer values, it turns out that a
      // *larger* SVD tolerance actually gives a better result than trying to pick a tolerance
      // close to machine precision.
      perf_log.push("svd_solve()");
      DenseVector<Number> ls_soln;
      A.svd_solve(rhs, ls_soln, /*rcond=*/1.e-10);
      perf_log.pop("svd_solve()");

      // Compute the difference between soln and ls_soln.
      perf_log.push("estimate error");
      DenseVector<Number> soln_error = soln;
      soln_error -= ls_soln;
      Real error_norm = soln_error.linfty_norm();
      perf_log.pop("estimate error");

      // We've seen dgelss fail to get the right solution at times,
      // but it's not clear why, and Lapack does not report any
      // errors.  Perhaps it is a tolerance/RCOND issue?
      if (error_norm > 1.e-2)
        {
          libMesh::err << "ls_soln = " << ls_soln << std::endl;
          libmesh_error_msg("libMesh svd_solve() got the wrong answer!");
        }

      libMesh::out << "libMesh: ||soln-ls_soln||_{infty}: " << error_norm << std::endl;
    }
  else
    {
      libMesh::out << "Skipping solve test, A has nrows < ncols." << std::endl;
    }
}



void eigen_solve_least_squares(const EigenMatrixType & A_Eigen)
{
  PerfLog perf_log("Eigen Solve Least-Squares");

  // 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");

  // The time for this is almost negligible compared to the time to form the SVD.
  perf_log.push("Eigen Least Squares");

  // Construct a rhs vector by multiplying the system matrix with a
  // vector of all 1's.  Then perform an SVD solve using this vector
  // as a rhs, and compare the result to the known solution.  Other
  // methods for solving a least-squares system of equations with
  // Eigen are discussed in
  // http://eigen.tuxfamily.org/dox-devel/group__LeastSquares.html.
  //
  // Note: if the original matrix has more rows than columns this makes sense
  // (more equations than unknowns).  If the matrix has more columns
  // than rows, there is most likely no solution.
  if (A_Eigen.rows() >= A_Eigen.cols())
    {
      EigenVectorType soln(svd.matrixU().cols());
      for (unsigned i=0; i<soln.size(); ++i)
        soln(i) = 1.;

      // Compute rhs
      EigenVectorType rhs = A_Eigen * soln;

      // Solve the system of Eqns A*x=rhs for x in the least-squares sense.
      EigenVectorType ls_soln = svd.solve(rhs);
      // libMesh::out << "ls_soln = " << ls_soln << std::endl;

      // Compute the difference between soln and ls_soln.
      Real soln_error = (soln-ls_soln).lpNorm<Eigen::Infinity>();
      libMesh::out << "Eigen: ||soln-ls_soln||_{infty}: " << soln_error << std::endl;
    }
  else
    {
      libMesh::out << "Skipping solve test, A has nrows < ncols." << std::endl;
    }

  perf_log.pop("Eigen Least Squares");
}