[WIP] Add Randomized SVD for Eigen which is a fast low-rank approximation and singular value decomposition for big matrices.

This is just a minimal example of implementation of Randomized SVD reference Kazuya`s code,some member functions of class RandomizedSVD has not finished yet since we are still not sure if this algorithm would be useful for Eigen users.We also test the performance with code belows:

#include<Eigen/SVD>
#include<chrono>
#include<iostream>

using Eigen::MatrixXd;
using Eigen::VectorXd;
using namespace std::chrono;

/*
  Computes spectral norm of error in reconstruction, from SVD matrices.

  Spectral norm = square root of maximum eigenvalue of matrix. Intuitively: the maximum 'scale', by which a matrix can 'stretch' a vector.
  Note: The definition of an eigenvalue is for square matrices. For non square matrices, we can define singular values: Definition: The singular values of a m×n matrix A are the positive square roots of the nonzero eigenvalues of the corresponding matrix A'A. The corresponding eigenvectors are called the singular vectors.
*/
double diff_spectral_norm(MatrixXd A, MatrixXd U, VectorXd s, MatrixXd V, int n_iter=20) {
  int nr = A.rows();

  VectorXd y = VectorXd::Random(nr);
  y.normalize();

  MatrixXd B = (A - U*s.asDiagonal()*V.transpose());

  // TODO make this more efficient (don't explicitly calculate B)
  if(B.rows() != B.cols())
     B = B*B.transpose();

 // Run n iterations of the power method
 // TODO implement and compare fbpca's method
  for(int i=0; i<n_iter; ++i) {
    y = B*y;
    y.normalize();
  }
  double eigval = abs((B*y).dot(y) / y.dot(y));
  if(eigval==0) return 0;

  return sqrt(eigval);
}

/*
  Test randomized SVD by computing decomposition of large (random) matrix
  Empirically, I found that increasing the rank of the randomized SVD doesn't decrease the reconstruction error enough to make it worth the increase computation. 
  You can check that the algorithm works by setting rank_rsvd=rank(M)
*/
int main(int argc, char* argv[]) {
  std::cout << "Making large (1000x500) random matrix" << std::endl;
  Eigen::MatrixXd M = Eigen::MatrixXd(1000, 500);
  srand((unsigned int) time(0));
  M.setRandom();

  // Full SVD
  std::cout << "Full SVD: ";
  auto start = steady_clock::now();

  Eigen::JacobiSVD<Eigen::MatrixXd> full_svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);

  auto now = steady_clock::now();
  long int elapsed = duration_cast<milliseconds>(now - start).count();
  std::cout << elapsed << " ms" << std::endl;

  // Randomized SVD
  int rank = 10;
  std::cout << "Randomized SVD with rank " << rank << ": ";
  start = steady_clock::now();

  Eigen::RandomizedSVD<Eigen::MatrixXd> rsvd(M,Eigen::ComputeThinU | Eigen::ComputeThinV, rank);

  now = steady_clock::now();
  elapsed = duration_cast<milliseconds>(now - start).count();
  std::cout << elapsed << " ms" << std::endl;

  std::cout << "Reconstruction error for full SVD (zero): " <<
    diff_spectral_norm(M, full_svd.matrixU(), full_svd.singularValues(), full_svd.matrixV()) << std::endl;
  std::cout << "Reconstruction error for rand SVD: " <<
    diff_spectral_norm(M, rsvd.matrixU(), rsvd.singularValues(), rsvd.matrixV()) << std::endl;

  return 0;
}

Get the output:

Making large (1000x500) random matrix
Full SVD: 6913 ms
Randomized SVD with rank 10: 23 ms
Reconstruction error for full SVD (zero): 0
Reconstruction error for rand SVD: 29.6648