Re: [eigen] RFC: BasicSVD class

[Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>]

On Fri, 7 Dec 2012, Gael Guennebaud wrote:

> Indeed, this simple algorithm seems pretty fast and it could even be even
> faster once we have fully optimized SelfAdjointEigenSolver. However
> I'm skeptical about its accuracy since it is squaring the entries that is
> typically something you want to avoid when using a SVD solver.

The algorithm indeed seems to struggle with ill-conditioned matrices, 
which have both small and large singular values. The test program below 
uses a 6-by-6 matrix with doubles and gives the following output:

-------------
Singular values (exact):
     1   0.01 0.0001  1e-06  1e-08  1e-10

Singular values (with JacobiSVD):
     1   0.01 0.0001  1e-06  1e-08  1e-10

Singular values (with BasicSVD):
          1        0.01      0.0001 9.99981e-07 1.05497e-08 1.43198e-09
-------------

Of course this does not mean that it's useless; there are probably 
circumstances in which BasicSVD is good enough as Mark indicates, and it 
would make a nice addition. The power method by itself would already be 
very nice.

Jitse


-----------------------------------

#include "BasicSVD.h"
#include <iostream>

using namespace Eigen;

int main()
{
  // The singular values
  const int size = 6;
  VectorXd singvals(size);
  singvals << 1, 1e-2, 1e-4, 1e-6, 1e-8, 1e-10;

  // Generate two random orthogonal matrices
  MatrixXd Uexact = MatrixXd::Random(size, size).householderQr().householderQ();
  MatrixXd Vexact = MatrixXd::Random(size, size).householderQr().householderQ();

  // Generate matrix whose SVD is to be computed
  MatrixXd A = Uexact * singvals.asDiagonal() * Vexact;
  std::cout << "Singular values (exact):\n";
  std::cout << singvals.transpose() << "\n\n";

  // Compute SVD using standard algorithm
  VectorXd singvalsJac = A.jacobiSvd().singularValues();
  std::cout << "Singular values (with JacobiSVD):\n";
  std::cout << singvalsJac.transpose() << "\n\n";

  // Compute SVD using new algorithm
  BasicSVD<MatrixXd> svdBasic(A, size);
  VectorXd singvalsBas = svdBasic.singularValues();
  std::cout << "Singular values (with BasicSVD):\n";
  std::cout << singvalsBas.transpose() << "\n";
}