Re: [eigen] RFC: BasicSVD class |

[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]

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